Primordial helium recombination I: feedback, line transfer, and continuum opacity 
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Precision measurements of the cosmic microwave background temperature anisotropy on scales 
t > 500 will be available in the near future. Successful interpretation of these data is dependent 
on a detailed understanding of the damping tail and cosmological recombination of both hydrogen 
and helium. This paper and two companion papers are devoted to a precise calculation of helium 
recombination. We discuss several aspects of the standard recombination picture, and then include 
feedback, radiative transfer in He I lines with partial redistribution, and continuum opacity from 
H I photoionization. In agreement with past calculations, we find that He II recombination proceeds 
in Saha equilibrium, whereas He I recombination is delayed relative to Saha due to the low rates 
connecting excited states of He I to the ground state. However, we find that at z < 2200 the 
continuum absorption by the rapidly increasing H I population becomes effective at destroying 
photons in the He I 2 1 P° — 1 X S line, causing He I recombination to finish around z ~ 1800, much 
earlier than previously estimated. 

PACS numbers: 98.70.Vc, 95.30.Jx 
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I. INTRODUCTION 

Cosmological recombination occurs when the photon 
gas in the early universe has cooled sufficiently for bound 
atoms to form. As the free electrons become locked in the 
ground states of these atoms, the opacity from Thomson 
scattering drops, and the signatures of thermal inhomo- 
geneities in the recombination plasma begin to stream 
freely across the universe. These signatures reach us to- 
day in the cosmic microwave background (CMB). The 
best limits on the CMB temperature anisotropy over the 
sky, down to ~ 0.4° are provided by the Wilkinson Mi- 
crowave Anisotropy Probe (WMAP) In conjunction 
with other surveys, these data tightly constrain cosmo- 
logical model parameters [l|, 0] • 

Most future and recent CMB experimental efforts aim 
to measure the polarization or t emp erature anisotropy 
@, 1, H I 0, 1 1 M EL IS on smaller 

scales than those measured by WMAP. Small-scale CMB 
temperature anisotropy measurements further constrain 
the matter and baryonic matter fractions, Sl m h 2 and 
Qbh 2 [l7| . the spectral index of the primordial scalar 
power spectrum, n s , and its possible running, a s . Mea- 
surements of n s will be integral to the viability or even- 
tual rejection of a wide variety of inflation models [l8l.[l9j. 

The most conspicuous features in the CMB tempera- 
ture anisotropy on small scales are acoustic oscillations 
H3, [HI, US H| and Silk damping [Mil!. Silk damping 
is determined directly by the free electron abundance, 
which is set by the recombination history. Here, pho- 
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tons in the ionized gas will diffuse over a characteristic 
scale k]j cx ^fn^ (where n e is the free electron number 
density and ko is the damping wavenumber), exponen- 
tially damping photon perturbations with wavenumbers 
k > krj. If the free electron density is ovcrprcdicted by 
recombination models, then kjj is also overpredicted and 
Boltzmann codes will predict too much power on small 
scales. Measurements of n s using CMB data including 
the damping tail region will then be biased downward. 
In light of the WMAP 3-year analysis it is crucial 
to understand how n s differs from 1. While the differ- 
ences among recently published recombination histories 
are too small to be important for the WMAP n 3 mea- 
surement, the corrections are expected to be significant 
(at the An s ~ 0.02 level) for Planck [2(|, and presum- 
ably also for high-£ experiments such as ACT || and SPT 
Q . Successful interpretations of data from the next gen- 
eration of small-scale CMB anisotropy experiments will 
depend on a solid understanding of recombination. 

Fundamentally, the problem in cosmological recombi- 
nation is to solve consistently for the evolution of the 
atomic level occupations and the radiation field (which 
has both a thermal piece, and a non-thermal piece from 
the radiation of the atoms themselves) in an expanding 
background. The highly-excited states are kept close to 
equilibrium by the high rates interconnecting them, and 
the rate of formation of the ground state in both he- 
lium and hydrogen is dominated by the occupation of 
the n = 2 states and the rates connecting the n = 2 
states to the ground state. In both helium and hydro- 
gen, the decay channels to the ground state through the 
allowed transitions are dramatically suppressed relative 
to the vacuum rates by the optical depth in the gas. In- 
deed, in both systems, the two-photon decay rate from 
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the n = 2 S state (singlet in the case of He i) is compa- 
rable to the rate in allowed decay channels. This is the 
so-called "n = 2 bottleneck" (27j |. 

The primordial recombination was first investigated 
theoretically in the 1960s [28|, [29[ using a simple "three- 
level atom" (TLA) approximation. The TLA tracks the 
abundance of ground state hydrogen atoms (H I Is), ex- 
cited hydrogen atoms (assumed to be in Boltzmann equi- 
librium), and free electrons. The TLA allows for recom- 
bination to and photoionization from excited H I levels, 
and allows excited atoms to decay to the ground state by 
Lya (2p — > Is) or two-photon (2s — > Is) emission. Subse- 
quently, a substantial literature developed, testing some 
of the assumptions of the TLA and extending it to the 
problem of helium recombination [3(| EE H3, S, [HI, [H] . 
Seager et al. [27[ provide the current benchmark precision 
recombination calculation by simulating 300 levels in H I, 
200 levels in He I, 100 levels in He II, interactions with the 
radiation field in the Sobolev approximation, basic hydro- 
gen chemistry, and matter temperature evolution. They 
found that the three level model [lH, H^] with a "fudge 
factor" inserted to speed up H I recombination was an 
accurate approximation to their full multilevel atom so- 
lution. Their recombination code, Recfast [36|, is pack- 
aged into most of the CMB anisotropy codes in common 
use, and underlies the cosmological constraints from the 
CMB, including those recently reported by WMAP. 

There are several reasons why it is now timely to 
revisit the cosmolo gica l recombination. Recent work 
[S3, H [H, S3, HE E3 has led to suspicion that some 
pieces of the recombination problem, such as matter tem- 
perature evolution [39l. |42| , two-photon transitions from 
high-lyin g st ates [381 ] , the effect of He I intcrcombination 
lines [33, [40] , departures of the I sublevels of hydrogen 
from their statistical population ratios (43[, and stimu- 
lated two-photon transitions [13, SU are not completely 
understood or were absent in past work. Two-photon 
transitions and intercombination rates speed He 1 recom- 
bination by facilitating the formation of the ground state, 
and as described in [38[ these could constitute a serious 
correction He 1 recombination, pushing He 1 much closer 
to equilibrium than in traditional models [27j | . These 
considerations led to significant changes in cosmological 
model parameter estimates [26| and warrant further con- 
sideration. (The allowed, intercombination, quadrupolc, 
and two-photon transitions up to n = 3 are shown in the 
Grotrian diagram, Fig. [T]) In addition, some of the is- 
sues considered during the 1990s, such as the effect of H 1 
continuum opacity on He I recombination, were not fully 
resolved; we will show here that this problem is much 
more complicated than previously believed. Finally, the 
imminent prospect of precision CMB data at high t from 
ACT, SPT, and Planck has "raised the bar" for theorists 
and demands a much better understanding of recombi- 
nation than was needed a decade ago. 

This is the first in a series of three papers devoted to 
the subject of helium recombination; the companion pa- 
pers (Hirata and Switzer |astro-ph / 0702 144) and Switzer 
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FIG. 1: Formation of neutral helium: a Grotrian diagram (up 
to n — 3) for He I. Singlet (5 = parahelium) and triplet 
(S — 1 orthohelium) levels and higher-order transitions give a 
rich system of low-lying transitions. Marked are two-photon 
transition (light dashed lines from 2 1 So and 3 1 5o), (allowed) 
electric dipole transitions (solid lines from 2 1 P° and 3 P°, like 
the Lyman series in H 1), the intercombination lines (heavy 
dashed lines from 2 3 Pf and 3 3 P°), and quadrupole transitions 
(from 3 D2)- The dipole transitions are treated in Sec. IV Dl 
using a Monte Carlo method with partial redistribution, for- 
bidden one-photon lines are treated in Sec. IIVI using an ana- 
lytic method for complete redistribution. (The energy levels 
are not drawn to scale.) 



and Hirata astro-ph/0702145) will be referred to as Paper 
II and Paper III. We have not yet completed the solution 
of hydrogen recombination because the radiative transfer 
is more complicated, the treatment of the high-n levels 
is more important, and the required accuracy is much 
greater. We do not believe our existing code is accu- 
rate enough for that application (see the discussion in 
Paper III) . This paper describes the standard multi- level 
atom code, as well as the most important improvements 
we have made, namely inclusion of feedback from spec- 
tral distortions, semiforbidden and forbidden lines, H 1 
bound-free opacity, and radiative transfer in the He 1 lines 
in the cases of both complete and partial redistribution. 
Paper II describes several effects that we find to have only 
a small influence on recombination: two-photon transi- 
tions in He I, interfering damping wings, and photons 
in the He 1 continuum. Paper III discusses the effect of 
the isotope shift between 3 He and 4 He resonance lines, 
Thomson scattering, rare processes, collisional process, 
and peculiar velocities. Although we find these effects to 
be small, it is necessary to investigate them in order to 
be sure of this. Paper III also contains a summary of the 
remaining uncertainty in helium recombination and its 
implications for the CMB power spectrum. 

One unresolved issue that was discussed extensively in 
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the 1990s is whether hydrogen photoionization can de- 
stroy photons in the He I n}P° — l x 5 resonances (in 
particular, 2 1 P° — 1 1 5) and the He I l 1 ^ continuum, 
pushing He I recombination closer to Saha equilibrium. 
This effect w as ig nored in the earliest work on helium re- 
combination [3OI ] . however Hu et al. [H| argued that the 
effect was strong enough to force helium into Saha equi- 
librium. The argument is essentially that the timescale 
for the absorption of these photons is much less than the 
Hubble time. On the other hand, Seager et al. [27| argued 
that hydrogen photoionization was negligible because the 
rate of He 1 line excitation was much faster than the rate 
of H 1 photoionization. We show here that both of these 
arguments are too simple: the effectiveness of H 1 opacity 
depends not only on the optical depth per Hubble time 
from H 1, but also on the width of the He 1 line and the 
redistribution of photons within the line. We have there- 
fore supplemented our level code with a radiative transfer 
analysis of the He 1 lines that takes account of H 1 con- 
tinuum opacity and emission/absorption of resonance ra- 
diation, including coherent scattering. The result is that 
H 1 opacity has little effect at z > 2200, because even 
though the H 1 optical depth per Hubble time can be 
3> 1, the H 1 optical depth within the width of any indi- 
vidual He 1 line is small, so the standard Sobolev escape 
rate is unmodified and recombination proceeds similarly 
to Seager et al. [27|. At lower redshifts the exponen- 
tially increasing H 1 abundance begins to absorb photons 
out of the He 1 2 1 P° — l 1 ^ line, rapidly accelerating He 1 
recombination and leading to its completion by z ~ 1800. 

Leung et al. [39| argued that the heating due to the 
13.6 eV energy release per H 1 recombination would in- 
crease the matter temperature and thus slightly delay 
recombination. However this energy release couples in- 
efficiently to the matter because it is delivered almost 
entirely to the photons, which are not absorbed by the 
matter [4^]. This energy actually goes into the spectral 
distortion to the CMB, not the matter. In order to prop- 
erly follow the matter temperature during recombination 
it is necessary to directly keep track of the changes in 
the total energy in translational degrees of freedom, and 
the number of such degrees of freedom (which decreases 
during recombination); the ratio of these quantities is 
&Blm/2. This was done by Seager et al. [27|, and we do 
not believe any modification of the basic methodology is 
necessary. 

Finally, the feedback of radiation field distortions from 
resonances onto lower- lying resonances is expected to 
slow recombination rates [27] . For example, escaped ra- 
diation from Ly/3 in H I could redshift onto the Lya tran- 
sition and excite atoms; analogous effects occur in He 1 
(and in He II, although this case is less important for 
CMB anisotropics). This is included iteratively in the 
level code developed here. 

In this paper, we will develop the effects through sev- 
eral sections. The multi-level atom code is described in 
Sec. [TTJ We introduce H 1 continuum opacity and its ef- 
fects on radiative transfer between He 1 lines in Sec. IIIII 



Sec. II VI discusses an analytic approach to transport with 
complete redistribution and continuous (H 1) opacity. 
Sec. IVl expands this treatment to include partial redistri- 
bution using a Monte Carlo simulation. We conclude in 
Sec. IVII In Appendix [Al we discuss the atomic data used 
here; Appendix [B] describes issues with transport under 
complete redistribution; and Appendix [C] describes our 
implementation of the Monte Carlo method for solving 
line profiles with coherent scattering. Appendix [D] exam- 
ines the limiting cases in which the Monte Carlo method 
reduces to the method of Sec. IIVI and to the widely used 
Sobolev escape probability method. Finally, Appendix [El 
describes the handling of photons above the He 1 pho- 
toionization threshold (24.6 eV) in our code. 



II. A NEW LEVEL CODE WITH RADIATIVE 
FEEDBACK 

The standard recombination scenario discussed in this 
paper is a homogeneous, interacting gas of protons, he- 
lium nuclei, electrons, and photons in an expanding back- 
ground. The background dynamics (i.e. the Fricdmann 
equation) are allowed to include other homogeneously 
distributed components (e.g. neutrinos, dark matter, and 
dark energy) but in our treatment we only consider their 
effect on H(z), neglecting any direct interaction with the 
baryons and photons. Several extensions to this stan- 
dard scenario have been proposed but will not be con- 
sidered here: energy injection from self-annihilating or 
decaying dark matter, primordial magnetic fields, small- 
scale inhomogeneities or peculiar velocities, and others 
HE H 53, M IS M M, El El, H| . The underlying 
physics in the standard problem is known, but the solu- 
tion is complicated by the interactions between radiation 
and atomic level occupations, and the variety of atomic 
rates. Many atomic rates involved are time-consuming 
to calculate, or not known well (for example, collisional 
and intercombination rates in He 1.) Collisional processes 
contribute to a lesser extent because of the high photon to 
baryon ratio during recombination. These are discussed 
in Paper III, where we argue that they are unimportant 
for helium recombination; in this paper we restrict our- 
selves to radiative processes. 

Recombination calculations must ultimately meet the 
practical constraint that they be fast enough (running in 
a few seconds) to act as inputs to CMB anisotropy codes. 
The general procedure for developing a recombination 
code that meets this constraint is to first "over-simulate" 
the recombination history, and prioritize the effects. One 
then develops a practical method that encapsulates the 
important physics. Additional effects can be included as 
parameterized corrections to the model. The TLA ap- 
proximation has been practical for CMB studies to date, 
and a constant correction factor is sufficient to bring it 
into agreement with more complete methods [27| . CMB 
researchers have been lucky that such a fast approxima- 
tion is possible, but in the post-WMAP era it is not clear 
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whether the TLA approximation will be sufficiently accu- 
rate. Tighter experimental constraints on the high CMB 
multipolcs in the Silk damping region will need to be 
matched by confidence that the underlying recombina- 
tion physics is well understood. Here, we use a multilevel 
atom code similar to Seager et al. [13] to assess the signif- 
icance of several new processes relative to a base model 
with one set of cosmological parameters. This level of 
detail is essential to developing an understanding of new 
physical effects. Unfortunately, it leads to a code that is 
much too slow for direct incorporation into Boltzmann 
codes. Developing and extension to the TLA approxima- 
tion (or similarly fast model) for a range of cosmological 
parameters will be the subject of later work. 

Many of the important details in the calculations are 
wrapped up in the notation, so we present the major 
symbols used here in Table [IJ 



A. Summary of the method 

For brevity, we will only highlight the physical argu- 
ments of the multilevel code, and differences from Seager 
et al (2?| ■ The multilevel atom code tracks levels up to a 
maximum principal quantum number n max . In this pa- 
per, we use a smaller model for H I with 245 levels (up 
to n max = 200), He I, 289 levels (n max = 100), and for 
He II, 145 levels (n max = 100), resolving I sublcvels up 
to n = 10, and including quadrupolc and intercombina- 
tion transitions in the He I rates. (A detailed discussion 
of atomic data can be found in Appendix [X] In Paper 
III, we investigate the effect of changing n max and find 
that these can be neglected here.) Unless stated other- 
wise, we assume a ACDM cosmology with £lb = 0.04592, 
Q m = 0.27, il r = 8.23 x 10 -5 , zero spatial curvature, 
masslcss neutrinos, a Hubble parameter of h — 0.71, and 
present-day radiation temperature T T {z = 0) = 2.728 K. 
The fiducial fractional helium abundance (i.e. ratio of 
helium to hydrogen nuclei) is /h c = 0.079. The Hubble 
rate in such a cosmology is 

h(z) = iW^A + n m (i + + fi r (l + z) 4 . (1) 

The number density of hydrogen nuclei is given by [27ll36| 



tih(z) = 1.123 x 10" 



1 + 3.9715/h, 



■(l + z) 3 cm- 3 , (2) 



where 3.9715 is the ratio of atomic masses of 4 He and 1 H. 
We will use the photon phase space density N(v) to track 
the radiation spectrum instead of the specific intensity 
J v , because the former is conserved along a trajectory in 
free space whereas J„ decreases as the universe expands. 
The relation between these is 



2hv 



■J„. 



(3) 



Photoionization and recombination contributions to 
atomic level population dynamics arc discussed in Seager 



et al. [27| . and we follow their treatment here. The rate 
of change of the average occupation of an atomic level 
is given by a series of bound-bound and bound-free rate 
equations. The photoionization rate from some level i is 
given by 



Pi = % 



(4) 



where N(i>) is the photon phase space density, Vth is 
the photoionization threshold frequency from level i, and 
a ic is the photoionization cross-section from that level, 
as a function of frequency (implicit in subsequent equa- 
tions). The photorccombination rate density for sponta- 
neous and stimulated processes is 



8tt ( n. 



n e n c 



LTE 



a, lc v 2 [l + N{v)]e- hv/k * T ™dv, 

(5) 



where T m is the matter temperature. The prefactor is 
the Saha ratio of the occupation of the level i (in local 
thermal equilibrium, LTE) to the free electron density 
times the continuum state density, 



LTE 



2Trm e kftT n 



3/2 



-£LeW*BT» (6) 



where i labels a bound state of a species, and c labels 
the continuum state (eg. for hydrogen this is H n) of 
the species; gt and g c are the state degeneracies. The 
bound-free rate for a level i is then 



dx i 
at 



(7) 



Except for several transitions in Hc I where transport is 
calculated separately to include new effects, single pho- 
ton bound-bound rates are calculated in the standard 
Sobolev approximation with complete redistribution, 

dx Q 

— = A u ^P s [x u (l - —xiN+], (8) 

at gi 

where A u ^i is the Einstein rate coefficient connecting 
an upper bound state u to a lower bound state I, and 
Af+ is the phase space density of radiation on the blue 
side of the line. In the "original" version of the code 
(i.e. before feedback is included) this is taken to be the 
Planck distribution, J\f + = Nm = {e hv l k ^ - 1) _1 . In 
the Sobolev approximation, the rates are modulated by 
the probability that a photon will escape from the reso- 
nance, allowing the average occupation state of the gas to 
change. The probability is associated with the Sobolev 
effective optical depth, 



TS 



so that 



SnH(z)^ ul 



ft = 



xi- 



9u 
9l 



Ts 



(9) 



(10) 



5 



TABLE I: Symbols used in this paper. Units of "1" mean the quantity is dimensionless. 

Symbol Units Description Equation 

a 1 Voigt unitless width parameter Eq. (|C5|) 

Ai^j s _1 Einstein spontaneous one-photon decay rate for i to j 

/ co h 1 fraction of line photon absorptions resulting in coherent scattering 

fi 1 fraction of line photon absorptions followed by transition to level i 

/i nc 1 fraction of line photon absorptions resulting in incoherent processes /i nc = 1 — / co h 

gi 1 degeneracy of level i, not including nuclear spin 

H s _1 Hubble rate 

K cm 3 s Peebles ^-factor [28] K = \f in j8nH 

n e cm -3 electron density 

rii cm~ 3 density of atoms in level i 

riH cm -3 total density of all hydrogen nuclei 

M 1 photon phase space density Eq. (j3j) 

A/"c 1 N in equilibrium with the continuum opacity Eq. (|41|l 

AL 1 A/" in equilibrium with the line opacity Eq. (1221) 

Al 1 modification of A/L used with coherent scattering Eq. (|74|l 

Api 1 photon phase space density for blackbody distribution Api = \j(e hv l k ^ T ^ — 1) 

A/± 1 photon phase space density on the blue (+) or red ( — ) side of a line Eq. (I25|l 

M 1 photon phase space density averaged over the profile Eq. (I50|l 

Pc 1 probability of photon emitted in line being absorbed by H I Eq. (|53|l 

P csc 1 escape probability from the line (equal to Ps in Sobolev approximation) 

Pmc 1 prob. of photon in MC being lost by H I absorption or redshifting 

Ps 1 Sobolev escape probability Eq. (I10|l 

Q ergs -1 heating per hydrogen nucleus per unit time 

Rij s _1 transition rate from level i to level j 

T m , T r K matter (T m ) or radiation (P r ) temperature 

x 1 frequency relative to line center in Doppler units x = (y — z^ii ne )/A^D 

x e 1 abundance of electrons relative to total hydrogen nuclei x e — n e /nn 

Xi 1 abundance of the state i relative to total hydrogen nuclei Xi = m /nn 

xt t 1 total number of free particles per hydrogen nucleus 

a; cm 3 s _1 recombination coefficient to leveli 

Pi s _1 photoionization rate from level i 

Ti s _1 width of level i 

rn nc s _1 Lorentz width of the line 
Fq ergcm _3 s _1 heating rate per hydrogen nucleus 

Aub Hz Doppler width of line 

Tj c Hz -1 differential optical depth from continuous opacity Eq. Q270 

Ai-,j s~ spontaneous two-photon decay rate from i to j 
Aq ergcm _3 s _1 cooling rate per hydrogen nucleus 

^linc v u \ Hz frequency of line center; v u \ for specific upper and lower levels 

i/th,i Hz photoionization threshold frequency from level i 

£ 1 rescaled photon phase space density Eq. ([76} 

n Hz -1 probability distribution of photon frequency in Monte Carlo 

<Ji C cm 2 photoionization cross section from level i 

or cm 2 Thomson cross section 

T co h 1 Sobolev optical depth from coherent scattering r co h = /cohTs 

ri nc 1 Sobolev optical depth from incoherent processes r mc = /inc^s 

tll 1 optical depth from continuous opacity between lines Eq. (|30|l 

ts 1 total Sobolev optical depth Eq. (|9]) 

<f) Hz -1 atomic line profile 

X 1 photon-atom scattering angle Eq. (IC1I) 
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B. Matter temperature 

The matter temperature, T m , departs from the radi- 
ation temperature due to adiabatic expansion, Comp- 
ton cooling, free-free, bound-free, and bound-bound pro- 
cesses 
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Here, we write the heat exchange terms to 
emphasize the individual processes, kinetic degrees of 
freedom in the gas, and we expand the description of 
free- free processes in [53] ■ Nothing is fundamentally new 
or unexpected in this treatment, and, indeed, we find 
that the departure of T m from T r is negligible during 
He I recombination. The rate of change of the matter 
temperature is related to the heating/cooling processes 
and adiabatic expansion in the gas through 



3fcB£tot 



(Qcs + Qs + Q 



bf 



ibb 



2FT m , (11) 



where the "es" subscript denotes Compton cooling and 
the others denote radiative atomic processes between 
bound and free levels. An individual heat exchange pro- 
cess Q has the general form 



Q 



Tq-Aq 3. 

TT^totKBJm, (12) 

«H 2 



where Tq and Aq are the heating and cooling rates (in 
ergcm~ 3 s _1 ), and itot accounts for processes that mod- 
ify the number of kinetic particles in the gas. The Comp- 
ton (electron scattering) cooling term is 



4 k B (T T - T m ) 



Qcs = -4x e ccr T a R T i : 



where we have used the radiation constant 



(13) 



7r /cj|/15c J ?i . The free- free contribution to T m is the in- 
tegral of the free-free opacity over the radiation temper- 
ature blackbody distribution (heating) minus the matter 
temperature blackbody distribution (cooling), 



Qff 



8tt 

TIrC 2 



is 2 [Af Pl (v, T T ) - Npi(v, T m )]hva s (v)dv, 

(14) 

where A/pi is the blackbody distribution, and the length 
absorption coefficient as{v) in units of cm^ 1 is J55| 



ot«{v) 



4e 6 



2tt 



3m e hc V ik B m l 

x(1 _ e -WfeT, 



rn— 1/2 -3 

■T m 1 v n e nmi 
l )3$- 



(15) 



The thermally averaged Gaunt factor gg is given by 
Sutherland [561 ]. In principle there is an additional cor- 
rection due to free-free radiation from electron-He II and 
electron-He in collisions, however these other species are 
an order of magnitude less abundant than H II, and devi- 
ations from T m = T r are negligible for helium recombina- 
tion, regardless. Therefore we have not included helium 
free-free radiation in our matter temperature evolution. 



The bound-free contribution to T m is from energy ex- 
changed in photorecombination and photoionization, and 
due to the change in Xtot! 



Q 



E 



n H 2 



-ain e x c k,K-Lr> 

n H 2 



(16) 



The energy exchanged per bound-free process includes 
the heating due to photoionization, 

r 87r f°° 

-=Xi—»l a ic v 2 N(v)h(v - vth,i)dv, (17) 

and the cooling due to recombination, 

LTE 



A 



Q.i 



7r / m 



Uh c \n e n c 

/>oo 

x / g ic v 2 \\ +N{v)]h{v- vth,i)dv. (18) 

We have computed these results using the blackbody 
function for M{v), since it is the thermal CMB photons 
that are responsible for the photoionizations and stimu- 
lated recombinations in the Balmer, Paschen, etc. con- 
tinua of H I, and the analogous continua of He I. The 
potential pitfall in this assumption is that some of the 
photons emitted in the He I 21.2 eV line (2 1 P°~1 1 S) 
photoionize H I directly from the Is level. This provides 
an additional source of heating during helium recombina- 
tion that is not included in our code. However a simple 
calculation shows that it is negligible. The maximum 
amount of thermal energy that can be injected by such 
photons is 21.2 — 13.6 = 7.6 eV per photon. Since a frac- 
tion fti e /xtot ~ 0.04 of the matter particles during this 
era are heliums (He I or He n), this implies an injected 
energy of ~ 0.3 eV per particle, even assuming every he- 
lium recombination injects the maximum amount of en- 
ergy. For comparison the total energy in translational 
modes of the matter is (3/2)^3^ = 0.6 eV per parti- 
cle. Even at the very end of helium recombination, the 
Compton cooling time is 



^Compto 



(3/2)x tot fc B (r m -T r ) 

Qcs 

3x t otm e c 



Sx e aTa^Tf 



7 x 10 b s: 



(19) 



in comparison, the timescale for helium to recombine is 
several times 10 12 s. The maximum fractional change in 
the matter temperature is then 0.3eV/0.6eV multiplied 
by the ratio of the Compton cooling to recombination 
time, which is of order 10~ 6 . This can be neglected. 

The contribution from bound-bound transitions (Qbt>) 
is negligible through a similar argument. Bound-bound 
processes do not change x to tj so the term Qbb is deter- 
mined by the average amount of energy delivered to the 
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gas by a photon that is injected into (or redshifts onto) 
a resonance line. In general, the net heating rate can 
only be calculated by a detailed radiative transfer anal- 
ysis [57] . However, as we saw in the previous paragraph, 
the injection of energy changes T m at the 10 -6 level if 
each He I recombination releases 7.6 eV of energy into 
the gas. Since there is a total of 24.6 eV available per 
He I recombination, T m would change by only several 
parts in 10 6 even if all of this energy were released into 
the gas by resonance line scattering. Therefore we can 
neglect the bound-bound contribution to T m . 

Time derivatives for the occupations and temperatures 
are converted to redshift derivatives according to 



(20) 



Equations dZJ) , (JHJ) an d (fTTj) give a set of stiff equations 
for the level occupations, as in Seager et al. (2?| that can 
be solved using the semi-implicit extrapolation (Badcr- 
Deuflhard) method [58j ]. 

Our matter temperature evolution is shown in Fig. [2] 
The important conclusion is that the fractional tempera- 
ture difference (T r — T m )/T T is negligible throughout he- 
lium recombination (z > 1600). This is in accord with 
the methods and results of Seager et al. [27|- At low 
redshifts when hydrogen rccombines the matter and ra- 
diation temperatures fall out of equilibrium, but that era 
is not the subject of this paper. 
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FIG. 2: The fractional difference in the matter temperature 
relative to the radiation temperature, 1 — T m /T r . Note the 
sign: matter is cooler than radiation at late times because of 
its different adiabatic index (5/3 versus 4/3). During He I 
recombination for 1600 < z < 3000, the fractional difference 
is < 10" 6 . 



the lower-lying state, and suppress the overall formation 
rate of the ground state. This radiation is manifest as a 
spectral distortion to the thermal spectrum of radiation 
incident on the blue side of the lines. (The calculation 
shares some similarity to [59[ where we showed that the 
relic recombination radiation is enough to ionize most of 
the neutral lithium population in the era following re- 
combination, except here the radiation is fed back onto 
recombination.) 

Feedback between levels can be calculated in a number 
of ways. The most direct way is to simulate a frequency- 
discrctized radiation field along with the atomic levels. 
In this method, the evolution of the radiation field bins 
and the evolution of the occupation states of the atoms 
are solved for simultaneously. This is numerically very 
difficult. Because of the huge range of rate scales in the 
system and the large number of radiation bins that it is 
necessary to track, this method quickly becomes difficult 
to manage and ill-conditioned. For this reason, we use 
an iterative method to include feedback in the level code. 
This is both practical to implement, and accurate. 

The level code calculates the radiation distortions gen- 
erated by transitions from excited states to the ground 
state in a first pass, for each redshift step, for each 
species. The first pass assumes a black-body radiation 
with local distortions in the Sobolev approximation. In 
a second pass, we transport the distortion generated by 
the (i + l)th transition to the ith transition to the ground 
state of the same species. That is, we only transport 
the radiation from the next higher-lying ground-excited 
transition, in the same species. (Inter-species feedback 
between He I and H I is discussed in depth in subsequent 
sections.) The distortion is recalculated for each level in 
the second pass, and these are transported to the lower 
levels and applied to a third pass of the level code. The 
iteration continues in this way. Because the iteration step 
accounts for much of the total feedback effect, subsequent 
iterations give progressively smaller corrections, and the 
procedure converges rapidly. The typical fractional con- 
tribution of the fifth iteration is |Ax e | as 2 x 10~ 4 - we 
stop there. 

For hydrogen, the distortion is determined by a simple 
vacuum transport equation for one resonance. In helium, 
the problem is complicated by continuous opacity from 
H I photoionization. Here, distortion photons are ab- 
sorbed by H I in flight, and may not redshift down to the 
lower level. This is treated in Sec. IIII B[ 

The radiation phase space density is convenient for 
cosmological calculations because it is conserved along 
a ray, yielding the bound-bound transport equation for 
complete redistribution (here for H i), 



C. Feedback 

Photons that escape a high-lying resonance will read- 
ily redshift to a lower-lying line owing to the (nearly) 
negligible optical depth between lines. This will excite 



ON 
dv 



T S (j>(v)[N{v) - A/l] 
x u gi 



TS<P(f) 



H{v)- 



XWu 



(21) 



where A/L is the radiation phase space density that would 
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be in equilibrium with the line, 
■A/l 



^ x u gi 

xi{g u /gi) 



(22) 



[This is obtained by setting the de-excitation rate 
A u ^i(l + A/") equal to the excitation rate g u A u ^iM / gi] 
Here we have used the approximation that the upper level 
is significantly less occupied than the lower level, which 
is valid when the lower level is the ground state since the 
energy of the first excitation in H I, He I, or He n is many 
times k^T v . The transport equation has the solution 



N{v) = ^-Cexp 
xig u 



-ts 



(23) 



where C is the constant of integration. This constant 
is obtained from the initial condition that the far blue 
side of the line has some phase space density A/+ , which 
implies 



C 



x u gi 
xig u 



■A/L 



(24) 



Thus the phase space density A/L on the red side of the 
line is 



A/L = A/L 



x u gi 
xig u 



-A/L (1-e-^). (25) 



Notice that as rs becomes very large, then the phase 
space density on the red side of the line approaches Al ■ 
Because M is conserved during the transport between H I 
and He n resonances (i.e. there is negligible continuum 
absorption or emission), the phase space density on the 
blue side of the H I Lyn resonance (i.e. ls—np) is simply 
the phase space density on the red side of the Ly(n + 1) 
resonance at an earlier time: 

■A/L (Lyn, z) = A/L \Ly(n+l),z f ] , 



l-(n + l)' 



l-rr 



-(1 + *)-!, (26) 



2 ] is the ratio of line fre- 



where [1 - (n + 1)~ 2 ]/[1 - 
quencies. A similar result holds for He n. Because of the 
existence of H I continuum opacity during He I recombi- 
nation, this result does not apply to He I; we will treat 
the He I problem in Sec. IIIII 



III. HYDROGEN CONTINUUM OPACITY AND 
HELIUM RECOMBINATION 

One of the major issues in recombination physics dur- 
ing the 1990s was whether helium recombines in Saha 
equilibrium, or is delayed. On the one hand, studies by 
Seager et al. [13] and Matsuda et al. [3(1 HH found that 
there is an "n = 2 bottleneck" in which He I recombines 
slowly because the two-photon 2 1 5 — l 1 ^ transition is 
slow and the 2 1 P° — l 1 ^ line is extremely optically thick; 



hence an excited helium atom has only a low probabil- 
ity of reaching the ground state and is most likely to 
be reionized by the CMB. These studies thus found a 
slower-than-Case B recombination for He I (where the 
v}P° — \ X S transitions are optically thick [60] and pro- 
cesses that depopulate a level through r^P — l^S are 
"blocked" by absorption of the same quanta, on average, 
thus greatly suppressing electron capture followed by di- 
rect cascade to the ground state). On the other hand, 
there is some neutral hydrogen present during the helium 
recombination era, and Hu et al. [44| argued that this can 
speed up helium recombination by absorbing resonance- 
line and continuum photons that would otherwise excite 
or ionize He I. There appears to be no satisfactory ex- 
planation in the literature for why the more recent works 
do not agree, but the difference in CMB spectra is sev- 
eral percent at high I J2?J so it is essential that the issue 
be resolved. This section, as well as Sees. II Vl andlVl are 
devoted to this issue. 

There are fundamentally three ways that H I continu- 
ous opacity could speed up He I recombination: 

1. Hydrogen can suppress feedback in the He I lines 
by absorbing He I line radiation before it redshifts 
down to the next line and excites a helium atom. 

2. If the H I opacity is very large, it could directly 
absorb He I resonance line photons, thus increas- 
ing the effective line escape probability above its 
Sobolcv value. 

3. Sometimes a helium atom recombines directly to 
the ground state. In the case of hydrogen, such 
recombinations are ineffective because the emitted 
photon immediately ionizes another atom. How- 
ever in the case of helium, it is possible to produce 
a neutral helium atom but have the emitted photon 
ionize an H I atom instead of He I. This results in 
a net recombination of helium. 



We treat mechanism $T] in Sees. IIII Al and IIIIBI The 
problem of absorption of He I line photons by H I (mech- 
anism ^2J| is more complicated. The physical picture is 
outlined in Sec. IIII CI and it is split into two cases. For 
the helium intercombination and quadrupolc lines, there 
is negligible coherent scattering within the line since the 
upper level has allowed decays (and allowed pathways to 
other states) whereas re-emission of the photon to the 
ground state of He I is semiforbidden or forbidden. This 
case is the simplest to consider and it is treated analyt- 
ically in Sec. IIVI The other case is that of the allowed 
He I n}P° — l l S lines, in which coherent scattering plays 
a key role alongside incoherent absorption/emission pro- 
cesses and H I opacity in determining the line profile 
and net decay rate. This situation is treated via Monte 
Carlo simulation in Sec. [V] Mechanism =$J3J produces 
a negligible effect and so we discuss it in Appendix [E] 
We will develop the distinction between coherent scat- 
tering and incoherent processes much more carefully in 
Sec. [mC] and Sec.lVA"! 
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A. Continuum opacity 

The photons emitted in resonant transitions in He I 
from excited states to the ground state have energies 
above the H I photoionization threshold. The opacity 
from photoionization influences transport both within 
and between He I lines. In this section, we describe how 
the continuum opacity is calculated, and Secs.lIVI andlVl 
describe details of transport subject to continuous opac- 
ity. Throughout, we use rj c to represent the continuum 
depth per unit frequency 



(It 
dv 



continuum 



nnxi s a c ic 
Rv ' 



(27) 



where a c \ is the photoionization cross-section of neutral 
hydrogen, and x\ s is the ground state occupation frac- 
tion (the excited states have much lower occupation num- 
bers and lower photoionization cross sections, so we ne- 
glect them). Stimulated recombination to H(ls) can be 
neglected. The continuum optical depth is also slowly 
varying as a function of frequency for the He I ground- 
resonance transitions because the energies are above the 
H I photoionization threshold and single-electron atoms 
such as H I posess no multiple-excitation resonances in 
their cross section. 

In standard recombination theory, the neutral hydro- 
gen population is well-described by the Saha distribution 
at early times (z > 1700): 



2Trm e k B T n 



3/2 



oXHl/feTn, 



(28) 



where xhii = 1 — ^hi ~ 1- One concern in using this 
equation to estimate the effect on helium recombina- 
tion is that radiation from helium recombination could 
knock the neutral hydrogen population out of equilib- 
rium. This would tend to lower the neutral hydrogen 
population, and decrease the continuum opacity. In the 
worst case, each He I recombination photoionizes a hy- 
drogen atom, giving the characteristic time scale for ion- 
ization, ionize ~ ^Hi/iHd- The Saha relaxation time 
for this perturbation to fall back to Saha equilibrium is 
well-approximated by the time scale for perturbations to 
decay in the TLA approximation (neglecting the H I re- 
combination rate at these high redshifts), 



^Saha 



2s^ls 



A, 



(Kn K xm) 



x ir e 



(29) 



where K = Xf jYa /8'KH is the Peebles A"- factor [28[ and 
(3h is the effective hydrogen photoionization coefficient 
from n = 2 in the three-level approximation. These are 
calculated, and shown for comparison in Fig. [3] It is ev- 
ident that a perturbation to the neutral hydrogen pop- 
ulation caused by He I photoionizing radiation quickly 
relaxes back to the Saha evolution. 



The unitless quantity rj c Ai>r) (where Avd is the the 
4 Hc Doppler width) during the He I recombination era is 
shown in Fig. |U 
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(d(x(Hel))/dt)/x(HI) ionization rate (Saha) 
(d(x(Hel))/dt)/x(HI) ionization rate (this paper) 
HI relaxation rate 
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FIG. 3: The ionization/relaxation timescales of H I during the 
period of He I recombination, assuming each He I recombina- 
tion generates a photon that photoionizes a hydrogen atom. 
Here we consider two He I recombination histories: one in 
equilibrium and the history derived here (Fig. I12[) . The "H I 
relaxation rate" is the inverse-timescale t jj*^ (Eq. I29|) for H I 
to return to Saha equilibrium if its abundance is perturbed. 
In either He I history, the ionizing radiation from He I is not 
sufficient to push H I evolution out of Saha equilibrium. 

The energy separations between transitions from ex- 
cited states to the ground state in He I are typically 
much greater than the optically thick line width. Thus, 
radiative transport in He I can be thought of as tak- 
ing place through two phases. In the first, continuum 
processes influence transport within a line. This modi- 
fies the transition rates, and sets the escape probability 
from a given resonance, which may exceed the Sobolev 
value. This is described in Sec. IIII CI In a second phase, 
continuum processes influence the transport of radiation 
between resonances, as described in the next section. 
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FIG. 4: The continuum optical depth d,T c /du — rj c times the 
Doppler width of He I 2 1 P° — l 1 ^ as a function of redshift. 
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Feedback from transport between He I lines and 
continuous opacity 



Feedback: contributions to HI, Hel and Hell 



Section III CI addressed the feedback of a radiation dis- 
tortion produced by a higher resonance on lower-lying 
resonances in H I and He II. For H I and He II, this trans- 
port is in free space in the approximation that the reso- 
nances are spaced more widely than their widths. (This 
approximation is not entirely valid for H I at late times, 
however this is not relevant to helium recombination and 
will be deferred to a later paper.) In He I, the picture is 
more complicated because transport is subject to opacity 
from the photoionization of neutral hydrogen. There is a 
sufficient neutral population that, when integrated over 
the photon's trajectory, the feedback between levels can 
be significantly suppressed. This is true especially near 
the end of the He I recombination and beginning of H I 
recombination. 

The algorithm presented earlier to include feedback it- 
cratively can be easily modified to include feedback sup- 
pression between lines: calculate the distortions for all 
resonances, and in the next iteration, multiply them by 
a suppression factor before applying them to the lower 
line. Changing variables to redshift, the total depth of 
the continuum between lines is 



dz 



H 



l + z 



(30) 



Let Mi be the radiation field on the blue side of the 
lower line assuming there is no line-line optical depth 
Tll- Then the nonthermal distortion produced by the 
higher state is Mi minus the Planck spectrum (since at 
times earlier than z ~ 1600 H I is in Saha equilibrium, 
to a good approximation), which is suppressed by the 
line-line depth. The final radiation field just above the 
frequency of the lower-lying line is then 



M+ = (Mi — Mpi) e~ TLL +M 



pi j 



(31) 



where the Planck spectrum is Mpi = 1 / \e hh '< /kBTv - 1). 

The effect of feedback is shown in Fig. [5l and the effect 
of including tll is shown in Fig. [BJ 
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FIG. 5: A comparison of the effect of feedback of a spec- 
tral distortion produced by higher-lying states on lower-lying 
states in the same species after several iterations, relative to 
the reference model, for He II, He I, and H I, where Ax is the 
abundance with feedback minus without feedback. Here, con- 
tinuous opacity from hydrogen photoionization between He I 
is included, (discussed in Fig. [6j). In all cases, feedback has 
the effect of retarding the formation of the neutral species. A 
larger number of iterations (~ 5) are needed for He I recom- 
bination than for H I or He II. Here, the uppermost line is 
the first iteration, moving down with further iterations and 
better convergence. 




1800 2000 2200 2400 2600 2800 3000 



C. Transport within He I lines and continuous 
opacity: physical argument 

Previous analyses [13, EH identify neutral hydrogen as 
a potential catalyst that could cause He I recombination 
to proceed closer to Saha equilibrium. This is because 
photons locked in the optically thick He I n 1 P°~l 1 S lines 
can ionize neutral hydrogen. This removes the photons 
and thus prevents them from re-exciting a helium atom: 

He^P ) -» Hc(l 1 S')+7, 

H(ls) + 7 -> H++e". (32) 

The purpose of this section is to investigate this process 
and related processes in the He I intercombination and 



FIG. 6: The effect of hydrogen continuum absorption on the 
feedback between transitions to the ground state in He I. 
Feedback slows He I recombination, but becomes increasingly 
less significant as more hydrogen recombines, increasing the 
bound-free opacity and absorbing the distortion. 



quadrupole lines. We will only discuss the physical sit- 
uation here, leaving the detailed calculation to Sec. IIVI 
(for intercombination and quadrupole lines) and [V] (for 
permitted lines). In those sections the objective will be 
to evaluate the increase in effective escape probability 
-Pesc (Eq. nU|) due to H I. Before we do so, however, we 
must discuss the issue of coherent scattering and its im- 
plications for what is meant by the "width" of a He I 
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line; we also provide a simple physical explanation for 
why H I opacity becomes important when it does (i.e. at 
z ~ 2200). 



1. Coherent vs. incoherent processes 

Optically thick lines in cosmological recombination are 
usually treated by the Sobolev method (see e.g. Sec. 
2.3.3 of Ref. [27| for a good discussion). This method 
describes decays from an excited level u of atoms to a 
lower level I without regard to how the atom reached the 
excited state. This a good approximation in the case of 
forbidden lines, and it is also good for very optically thick 
lines in the absence of other sources of opacity (in this 
case the radiation field in the optically thick part of the 
line comes to equilibrium with the population ratio of u 
and I, and thereafter the details of the radiative transfer 
matter little). However in the case of He I recombination, 
where He I and H I compete for photons, no equilibrium 
is established and the resonance line profile matters. In 
particular, the frequency distribution of photons emitted 
in e.g. 2 1 P° — need not be a Voigt profile and indeed 
need not be the same as the absorption profile <f>(i>). The 
reason is that sometimes level u is populated by resonant 
absorption of a photon from level I, which subsequently 
decays back to I without any intermediate interactions. 
In this case, conservation of energy requires the frequency 
of the final photon to be the same as the frequency of the 
initial photon in the atom's rest frame, and hence such 
events will be called "coherent scattering." In the co- 
moving frame the photon's frequency undergoes a small 
fractional shift of order v/c, where v is the atomic ve- 
locity. This small frequency shift has a minor influence 
on recombination, and will be included in the calculation 
of Sec. El For emphasis we will continue to call these 
scattering events coherent with the understanding that 
coherence is meant to be exact in the atom's frame; this 
should not lead to any confusion since there is no type of 
scattering that is coherent in the comoving frame instead 
of the scatterer's frame. (Note, however, that the term 
"coherent scattering" applied in the comoving frame does 
appear in literature.) 

The process described above differs from "incoherent 
scattering" [6l|, [62j , in which the excited atom undergoes 
other interactions (almost always involving one or more 
photons) before returning to I. It is only in the latter case 
that the final photon's frequency distribution can be de- 
scribed by a Voigt profile (complete redistribution) . It is 
also possible for the atom in the excited state u to be- 
come ionized, which we will consider to be an incoherent 
process since if the electron later recombines (probably 
onto another atom) the spectrum of emitted radiation 
will bear no memory of the frequency of the photon that 
initially excited the atom. Based on this distinction, one 
may split the Sobolev optical depth for a line into pieces: 
ts = T coh + rinc, depending on the fractions / coh and / inc 
of photon absorptions in the line that go to coherent or 



incoherent processes, respectively. Note that by classi- 
fying ionization from the excited state as an incoherent 
process, we ensure f mc + / co h = 1. 

The practical implication of this distinction is that 
when solving for the phase space density across 
an optically thick line, incoherent processes can play a 
dominant role even if f- lnc <§C 1 , as is the case for the He I 
2 1 P°-\ 1 S line. This is because incoherent scattering can 
transport a photon from the line center to a far damping 
wing (or vice versa) in one scattering event, and an exci- 
tation followed by ionization can remove a photon from 
the resonance line. In contrast, coherent scattering can 
only move a photon far into the damping wings by tak- 
ing many "baby steps" using the Doppler shift from the 
atom's velocity, and it by itself cannot create or destroy a 
line photon. We distinguish two physically distinct cases 
here: one where the coherent scattering is negligible (to 
be studied in detail in Sec. IIV|) and a more complicated 
case of partial redistribution where it is not (Sec. W\ . 

2. When is continuum opacity important? 

The He I 2 1 P° — l^S line is in the extreme UV at 21.2 
eV, where the optical depth from neutral hydrogen pho- 
toionization is nearly constant in the neighborhood of the 
resonance. The continuous opacity from neutral hydro- 
gen is also present in radiation transport within the He I 
n 1 P°-l 1 S series, [He i] r^D-V-S, and He i] n i P -l 1 S. 
Seager et al. [27j concluded that this effect is negligible for 
2 1 P° — l 1 ^ because the neutral hydrogen photoionization 
rate from photons in the 2 1 P° — l 1 ^ resonance is orders 
of magnitude lower than 2 1 P° — 1 : 5 photoexcitation rate 
during the He I recombination history, due to the sparse 
population of neutral hydrogen. A better criterion for the 
significance of neutral hydrogen is whether or not contin- 
uous opacity affects radiative transport within the line. 
The two natural scales in the transport problem are: (i) 
the frequency 77" 1 that a photon can be expected to tra- 
verse by redshifting before it is absorbed in a hydrogen 
photoionization event, and (ii) the range of frequencies 
A^ii n o over which the line is thick to incoherent scat- 
tering/absorption. If rj~ 1 3> Ai^inc then helium atoms 
will re-absorb the resonance radiation from other helium 
atoms, without regard to the H 1 abundance. If, on the 
other hand, r/" 1 < A^i; no , then H 1 can destroy pho- 
tons that would otherwise have re-excited helium atoms, 
and thereby accelerate He 1 recombination. We thus care 
about the continuum optical depth within a line, 

r c = r? c A^i ne . (33) 

In the case of lines that are optically thick into the damp- 
ing wings, such as the He 1 n x P° — \ X S lines, Afu ne may 
be calculated by integrating the asymptotic line profile 
<j){v) ~ T\i nc /A-K 2 /S.v 2 until the optical depth in the wing 
becomes unity: 

A Inline 7"inc . . 
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Here rii no is the Lorentz width of the line, Ti nc is the 
Sobolev optical depth through the line from incoherent 
processes (i.e. all absorption processes other than coher- 
ent scattering; this will be calculated in Sec.|V|, and r) c is 
the differential continuum optical depth, Eq. (|27|) . Fig. [7] 
compares the optically thick linewidth (from incoherent 
processes) to the inverse differential optical depth rj~ l 
and the Doppler width for the 2 1 P° — l x 5 transition in 
4 Hc. The line may be optically thick in the Sobolev sense 
out to much larger frequency separations due to coher- 
ent scattering, but a coherent scattering event results in 
no net change in the atomic level populations and hence 
does not directly affect recombination. (It only has an 
indirect effect by changing the radiation spectrum.) 



Frequency scales in the 2 1 P°-1 1 S line 




1600 1800 2000 2200 2400 2600 2800 3000 



FIG. 7: Inverse of the differential optical depth rj c from hy- 
drogen photoionization as a function of redshift, compared to 
the optically thick line-width due to incoherent processes in 
the 2 1 P° - l 1 ^ line. Continuum processes start to become 
important over scales inside the (incoherent) optically thick 
part within the line around z — 2100. Also plotted is the 
Doppler width of the line, emphasizing that the line is opti- 
cally thick out into the wings. Continuum processes do not 
act on scales smaller than the Doppler core until z < 1800. 

The general problem of the escape probability includ- 
ing the continuum opacity and coherent scattering as 
well as incoherent emission/absorption processes is quite 
complicated. Therefore we will solve it in two steps. In 
Sec. IIVI we will solve the problem without coherent scat- 
tering. This is a conceptually simple problem - all one 
has to do is compute the probability of a photon either 
redshifting out of the line or being absorbed by H I be- 
fore being re-absorbed by He I - and the machinery for 
solving it has already been developed for the theory of 
line transfer in stellar winds [63j | . Despite its simplicity, 
the solution in Sec. IIVI is an accurate description of the 
intercombination and quadrupole lines because there is 
negligible coherent scattering in these lines. (The reason 
is that if an atom absorbs a photon in these lines and 
reaches a n 3 P° or n Y D level, its next step is almost al- 
ways to undergo one of the allowed decays rather than to 



emit a photon in the intercombination or quadrupole line 
and go back to the ground state.) In Sec. [V] we address 
the problem with coherent scattering as well as incoher- 
ent processes and continuum opacity, and describe its 
solution via a Monte Carlo simulation. 



IV. THE MODIFIED ESCAPE PROBABILITY 
WITHOUT COHERENT SCATTERING 

This section is concerned with calculating the net de- 
cay rate in He I lines with continuous opacity (from H i) 
but no coherent scattering. This is a line radiative trans- 
fer problem with "complete redistribution" in the sense 
that a photon emitted in the line has a frequency distri- 
bution given by the intrinsic line profile, independent of 
the spectrum of incident radiation: 

P(foutl^in) = <H^out)- (35) 

There is negligible coherent scattering in the He i] 
n 3 P o _ 1 i 5 and j He jj n i D _ x i s lines ( see S e C . |rnC]) . 

and results in this section are readily applied to those 
lines. The macroscopic picture is that in emission, a 
line in the gas always has the same shape, regardless of 
the excitation field. Throughout, is Voigt-distributed 
and accounts for the natural resonance width of the line 
and Doppler broadening from thermal motion, set by the 
matter temperature. In Sec. [V] we discuss the case of a 
line in a gas where some fraction of the photons are re- 
emitted coherently (with the same frequency as the in- 
coming photon in the rest frame of the atom) and some 
arc rccmittcd incoherently. 

The fundamental quantity that determines the transi- 
tion rate between a lower level / and upper level u is the 
radiation phase space density integrated over the atomic 
absorption profile, N(v u i)- In the case of complete redis- 
tribution through incoherent scattering in the presence 
of continuous opacity, the radiation phase space density 
evolves as 

M = A/Hubble + Af cont + -Wine. (36) 

(We assume that the line width is small compared to v u i 
so that factors of v /v u i appearing in the transport equa- 
tion can be dropped.) In this section, we will develop 
each of these terms and solve the transport equations for 
the radiation profile over the line and the escape proba- 
bility. Of these terms, the Hubble redshifting term is, 

-^Hubble = Hv u i——, (37) 
av 

We will work in the steady-state approximation, where 
M = 0; possible corrections to this are considered in Pa- 
per II. Then, frequency provides a convenient domain 
over which to solve for Af, moving the Hubble term and 
dividing by Hv u i . Breaking up N CO nt and Mnc into emis- 
sion and absorption pieces, there are four terms that 
appear in the steady-state equation for the phase space 
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density: absorption and emission by continuum processes 
and emission and absorption by incoherent processes in 
the line. We may write these as 



dN 
dv 



dN 
~dv~ 

( dN 
V dv 



cont — abs 



dN 

( dN 
V dv 



inc — abs 



(38) 



The continuum absorption term has already been de- 
termined (Eq. I27|) and is rj c N(v). The line absorption 
term depends on the line profile and is Tg(p(v)N(v) since 
ts4>(v) is the optical depth per unit frequency. The con- 
tinuum emission term is 



\ / cont — em 



-VcNc, 



(39) 



where Nc is the phase space density of photons that 
would be required for the reaction 



H(ls)+7 



H H 



(40) 



to be in equilibrium. The usual equilibrium constant 
argument shows that Nc oc nHii^e/^m; the proportion- 
ality constant can be determined from the principle of 
detailed balance as 



N c (v) = 



n e n c 

Tit 



LTE 



(41) 



Since we have seen that hydrogen is in Saha equilibrium 
during the relevant era, and that T m = T r , we have 
N c {v) = e - hv/k * T *. The He i] n 3 P° - l 1 ^ and [He i] 
n^D — l^S lines are at very high energies (> 20.6 eV), so 
the photon phase space density < 1 and we can ne- 
glect stimulated emission and other consequences of the 
photon's bosonic nature. The line emission term is 



( dN 
V dv 



= -t s 4>(v)Nl, 



(42) 



where A/l is the phase space density that would ex- 
ist if only the line processes were important. This 
can also be determined by setting the excitation 
rate xi(g u /gi)A u ^iNh equal to the de-excitation rate 
x u A u ^i(l + JV L ): 



N L = 



xi(g u /gi) - xu 

The transport equation is thus 
dN 



(43) 



Vc (N-N c )+t s cP(v)(N~N l ). (44) 



It has the general solution: 



N(v) = e* l(iy) \ C - I [Ncv.+N^Ts^e-^dvj, 

(45) 



where 



$i(f) = / [r? c + Ts<p{v)]dv. 



(46) 



Here v\ is an arbitrary but fixed frequency, and vi is 
the frequency at which we set the initial condition. It is 
convenient to expand Eq. (|45|) into the pieces that depend 
linearly on the constant C and on Nc and Nl: 

N(v) = Ch{v) +N c I C (v) +N L I L (v), (47) 

where the individual profiles are 

h{v) = e* 1 ^ 

/ ffce-* 1 and 



Il{v) 



{v)e-* l[0) dv. (48) 



We integrate this phase space density over the profile, 
and break the integral into three pieces that emphasize 
the physical processes: 



N = Ch + N c ic + NJ L - 



(49) 



Here the overbar denotes averaging over the line profile, 
e.g. 



N = / <t>{v)N{v)dv. 



(50) 



Bringing the outer exponent under the integral in the 
expression of II'- 



II = - 



x cxp ■ 



V) I T S (j)(v) 



[Vc + Ts<f>{y)]dy > dvdv. (51) 



We will take the starting frequency to be on the far blue 
side of the line (z^ — so), appropriate for expanding me- 
dia [Hj]. With this choice, it is convenient to switch the 
order of integration and absorb the leading minus sign. 
The II line integral is related to the original Sobolcv 
problem by setting rj c = 0, 



1 



1 



1-Ps 



Vc=0 



TS 



(52) 



Simply, Ps = 1 — Il(Vc = 0). We follow Rybicki and 
Hummer [63j in then calculating the difference in the line 
integral with and without continuum absorption (which 
is related to the probability of absorption between inco- 
herent scattering events, Pc)' 



AI L = It. - It 











/ * 






J — oc 


J Vl 



x exp I -ts / (p{y)dy 



- 1 



dvdv. 

(53) 
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We integrate this numerically over a Voigt profile to find 
the probability of absorption by a continuum process. 
The Voigt profile and its integral have a wide dynamic 
range-details of the evaluating the integral are described 
in Appendix [B] The main goal is to find the integral of 
the radiation phase space density over the line profile (AT 
in Eq. I49p , and the modified escape probability. Exam- 
ples of line radiation profiles with incoherent scattering 
and continuum opacity are shown in Fig. [8] 

The overall value of N is required in order to compute 
excitation and de-excitation rates. With the boundary 
i>2 — * oo, we have $1(^2) - * 00 and hence C — ► 0, while 
Ii remains finite at fixed v\. Therefore 



from which we find 



N = N C I C + MJ L - 



(54) 



Now, if we had Ac = A"l, the solution to Eq. l|4"I)) would 
be simply Af = A/"l and hence A7 = A/"l- Thus Tc+h = 1 
and hence 



N = Nclc + A/l(1 - Ps + A/i) 
= Af c (P S - AI L )+Af h (l-P s 



AJ L ). (55) 
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FIG. 8: The radiation profile near the intercombination 
line I 1 So — > 2 3 P° resonance with Voigt parameter a — 
riinc/(47rA^D) = 10~ 5 for the Sobolev optical depth rs = 2.8 
and Nc = 0, for several sample continuum optical depths 
?7cA^d, showing the effect of continuous opacity. Because of 
the low optical depth and small natural line width, very little 
radiation extends more than three doppler widths above the 
line, and the effect of the continuum is significant in relaxing 
the radiation phase space density near line-center. 



The net downward transition rate from u to I is 

flu 



$ I line — A u - 



9i 



(56) 



For x u <§C xi and Af « 1, this can be re-expressed in 
terms of AL as 



£|line = A u ^i—xi(ATl - A"), 
9i 



(57) 



illtoe = A u ^x,(P s - AJ l )(Al - Ab) 
9i 



9u 
9l 



xiNq 



(58) 



where P csc = Ps — A/^. Note that AIl is manifestly 
negative according to Eq. (|53[) . so continuum opacity en- 
hances the escape probability. The transition rate from 
from level u to I is then set by the occupations of the 
states, and by the radiative escape probability, P esc . 

Before continuing, we note one subtlety of our analysis. 
We have assumed that the continuum on the blue side of 
the line is optically thick to H 1 photoionization, or equiv- 
alently, we neglected feedback from unabsorbed spectral 
distortions from a higher-frequency line. [Formally this 
was done when we argued that $1(^2) — > 00 ■] This is not 
as restrictive an assumption as it seems (even though it is 
violated at the beginning of He 1 recombination), because 
the widths of the He 1 lines are small compared to their 
separation. Thus if the continuum opacity between lines 
is not large, then the continuum opacity within a line 
can be neglected, in which case we know that the correct 
escape probability is the Sobolev result, P osc = Ps (and 
the converse, where opacity is important within the line 
only once with transport of the distortion between lines 
is suppressed) . Since our solution here for P esc reduces to 
the Sobolev optical depth in the limit where continuum 
opacity within the line is turned off, it follows that using 



i|linc — Pcsc-A u ^l ( X u X[Af+ 

91 



(59) 



i.e. replacing Nc in Eq. (jS"8")) with Af+, will recover 
the correct solution in both the case where continuum 
opacity within the line is important (where we have 
A/+ = Afc) and the case where it is not. 

The total escape probability is too time consuming to 
calculate on a case-by-case basis, as part of the level code. 
Instead, we calculate a lookup tabic of modified escape 
probabilities as a function of the He 1 ground state frac- 
tion and rcdshift, assuming that the neutral hydrogen 
population is in Saha equilibrium during He 1 recom- 
bination. The details of this calculation are given in 
Appendix [Bj The Monte Carlo methods developed in 
Sec. IV Dl also generate a grid of modified escape prob- 
abilities as a function of redshift and the occupation of 
the He 1 ground state. These arc log-interpolated and 
applied in to the recombination level code, as described 
in Sec. Ell 



V. THE MODIFIED ESCAPE PROBABILITY 
INCLUDING COHERENT SCATTERING 

In the previous section, we considered the problem of 
calculating the decay rate of He 1 lines with no coherent 
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scattering. This section aims to include coherent scatter- 
ing in the problem, which is necessary in order to handle 
the Lyman-like series of lines, He I n 1 P° — l 1 S. Coherent 
scattering complicates the problem because the frequency 
distribution of photons emitted in the line depends on 
the existing spectrum of radiation (photons are not com- 
pletely redistributed across the line). Our approach to 
this problem will be to realize that a helium atom reaches 
the upper level n 1 P° by emitting or absorbing a photon, 
and then it leaves this level by emitting or absorbing a 
photon. Therefore, all processes involving this level can 
be regarded as a form of resonant two-photon absorption, 
resonant two-photon emission, resonant Raman scatter- 
ing, and resonant Rayleigh scattering. Physically one 
should write down the rates for these processes and solve 
the relevant level population/radiative transfer problem. 



A. Setup 



level.) The overall width of the level u is then 

b>u 9u 



(60) 



(Note that the summation over b includes an implied inte- 
gration over continuum states.) It is convenient to define 
the rates 

Rui = f +N{v ui )\ Ei < E u 

\ Ai^ u {gi/g u )M{v iu ) Ei > E u 

so that T u = A u ^i + Rui,i^i- The fractional con- 
tribution to this width from a final level will be denoted 
fi = R U i/T Ul where we identify R u i = A u ->i. This should 
be interpreted as the probability that u will transit to i, 
given all of its options. The radiative rate for one-photon 
radiative transitions from some level i ^ I to u is 



The goal here is to replace the usual treatment of the 
line through one-photon processes by an inherently two- 
photon treatment (two-photon absorption, two-photon 
emission, and scattering) in which u represents the inter- 
mediate state. In the case of each He I n 1 P° — l 1 S line, 
we will denote the lower level by I = l x 5 and the upper 
level by u = n l P°. The rates for two-photon processes 
in the vicinity of resonance can be expressed in terms of 
the constituent one-photon absorption and emission pro- 
cesses, and depend in particular on the branching frac- 
tions that determine the fate of a helium atom in level 
it. Throughout, we will make several assumptions about 
the radiation field and rates during recombination: 

1. Af(Uyx) <C 1 implies that stimulated emission in 
He I n 1 P°-l 1 S lines can be neglected. (This is 
because hv u i ksT r so the Wien curve has Af <C 
1, and the spectral distortion raises Af to at most 

~ Xnipo/Zx^s < 1.) 

2. The transitions from u to other excited states 
(or the continuum) see approximately a blackbody 
spectrum and are Sobolev optically thin. 

3. The time for x u to change significantly is much 
longer than the lifetime of the state, so we can work 
in a "steady state" . 

An atom in level u has three possible fates. It could 
decay to the ground level I with the rate A u —,i (we ne- 
glect stimulated emission in the 2 1 P° — l 1 ^ line under 
assumption 1). It may also decay to another level a with 
lower energy E a < E u , with rate A u ^ a [l + Af{v ua )]. A 
third possibility is that an atom in level u could absorb 
a photon and transit to an even higher level b with rate 
Ab—> u (gb/ 9u)N(i>bu)- (Note that b could be a continuum 



Rim = 



Au-tiQju/ 9i)N(v u i) Ei < E. u 
Ai^ u [l+N{v lu )\ Ei>E v 



(62) 



The radiative rate for one-photon excitation from I to u 
is similar, except that one must average the photon phase 
space density over the line profile because the phase space 
density may be strongly frequency-dependent (unlike the 
u «-> e transitions where one is dealing with a blackbody 
radiation field): 



R lu = A u ^M. 

9i 



(63) 



Figure [9] summarizes the rates P; u , Ri u , Ri u , and Ri u for 
the lower excitation states of helium. 

It is straightforward to write down the net rate of pro- 
duction (or destruction) of each level via two-photon pro- 
cesses involving the excited state u, in terms of the rates 
(Eas. 1501 - 153"]) . It is simply the difference of destruction 
rate of level i via transition to it, and the rate of produc- 
tion of i from all modes involving u: 



Riu^i ~l~ ^ ^ •KjRjufi 



Riu^i ~l~ Rui^u ^ ] XjRju- (6^) 

If we knew exactly the rate coefficients R u i and Ri u , 
Eq. ([64]) would be easy to incorporate into the level codes. 
That is because we can simply write the steady-state pop- 
ulation of u as its production rate multiplied by its life- 
time r^ 1 , 



z u — T u XjRj U . 



(65) 



Then for i ^ l,u, Eq. (J5J} is exactly equivalent to the 
standard rate equation (Eq. |5J). 
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FIG. 9: The transition 2 1 P°-1 1 S interpreted as a two-photon 
process with 2 1 P° as an intermediate resonant state. In an 
incoherent scattering through He I l^So <-> 2 Pf, the incom- 
ing photon excites l 1 So — » 2 1 P° . The atom absorbs a second 
photon and explores several intermediate states before decay- 
ing through 2 1 P° - l 1 ^, emitting a 2 1 P° - l 1 ^ photon with 
complete redistribution. In coherent scattering, only l 1 ^ 
and 2 1 P° are be involved, and the outgoing photon is emit- 
ted with the same energy as the incoming photon in the rest 
frame of the atom. 



The rate equation for the ground level is modified to 
become 



xi\ u = -RluXl + R u ix v 
and that for u becomes 

= — RidXu + Ri 



(66) 



(67) 



These two are also similar to Eq. ([5]) except that the left- 
hand side of Eq. (|67j) is zero, as appropriate for a steady- 
state solution, and the rates R u i and Ri u may differ from 
the Sobolev values. Indeed, the only modification needed 
in the level code is that these rates need to be replaced 
with the values A u ^i and Eq. (|63p . The only nontrivial 
part of the calculation is to compute Af appearing in 
Eq. (|63|) . The calculation of Af will occupy the rest of 
this section. 

In some parts of this section we will define the addi- 
tional branching fractions / co h = /; and /j nc = 1 — 
These represent the branching fractions for absorption of 
a line photon to result in coherent scattering (/ co h) or in- 
coherent scattering (/mc)- We also use the portions of the 
Sobolev optical depth due to these processes, r co h = T sfi 
and r inc = rs(l - fi). 



B. The equation of radiative transfer 

In order to compute Af, we need to construct and solve 
the radiative transfer equation for the photon phase space 
density in the vicinity of the line, Af(v). The line pro- 
file evolves according to four effects: the Hubble red- 
shifting of the photons; absorption and emission in the 



H I continuum; coherent scattering; and incoherent emis- 
sion/absorption processes (i.e. resonant Raman scat- 
tering and two-photon emission/absorption, whose line 
profiles do not depend on the incident radiation field). 
Schematically, we write 

Af = A/Hubble + Kont + Koh + Wine • (68) 

The continuum contribution was solved in Sec. IIVI 



Af cont = -Hu u ii] c {Af -Af c )- 



(69) 



A photon has a probability per unit time of undergo- 
ing a coherent scatter given by Hv u iT co h<f>(v), where 
T C oh = Tsfi, the Sobolev optical depth to any type of 
absorption times where /; is the fraction of photon 
absorptions that are followed directly by emission, <j){v) 
is the line profile in units of fraction of the integrated 
profile traversed per unit time (s _1 ), and Hv u i converts 
this to the fraction of the integrated profile traversed per 
unit frequency (Hz" 1 ). Thus we have 

Af co h(v) = -Hv ul T coh (f>(is)Af{is) 

+Hu ulTcoh J <f>{v')N{v') V {v\v')dv' . (70) 

Here p(v\v') is the probability distribution for the out- 
going frequency v of a coherently scattered photon con- 
ditioned on its ingoing frequency v' . (This is commonly 
known as a redistribution function, and the relevant case 
here is of the R\\ type [65|.) It satisfies J p(v\v') dv = 1. 

Finally we come to incoherent processes, A/inc- The 
probability per unit time of incoherent scattering (i.e. 
excitation of an atom to u followed by transit to a state 
other than I) is Hv u iT- mc <j){v). The rate of incoherent 
emission processes (two-photon emission or resonant Ra- 
man scattering, i — > u — > I, with i ^ I) per H nucleus per 
unit time is a^-Rm/z- Thus we have 



Mnc(^) = -Hu u iT- inc (j>{v)N{v) 



nnc 
~8tw 2 , 



(71) 



(Here 7iHC 3 /87ri/^ is a conversion factor: it is the number 
of H nuclei in a volume containing one photon mode per 
unit frequency at v u i-) We note that /; = A u ^>i/T u , and 
using x u <C xi , we may write 



fl 



Hv u iT S 



Hv u igiT S 



xi(g u /gi) - x u 



r- 1 



(72) 



This turns Eq. ([TTJ) into 

M nc M = Hv ul T inc cj>(v)[Afl 0) - (73) 
where we have defined 



Af?* = 



9l 



XiRiu • 



(74) 
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Note that Af^ is not necessarily equal to the phase space 
density of photons in equilibrium with the line, which 
would be Al = giX u /g u Xi for x u <C xi. However, the two 
are related. 

In steady state one may combine Eqs. (j3"7| . (|6"9")l . ([70)) 
and (J73J) and set M -> to yield 

^ = r^A^-A/c] 



, ^v)M{v) - J <f>(v')M(v')p(v\v') dv' 

+t^{v)\M{v)-M { *\ (75) 

The parameters Mc and A/^ ' may be eliminated from 
this equation via the introduction of the variable 



A/» - A/p 
M^-Mc 



A^)=A/c + (A^ 0) -A/c)^). 
Expressing Eq. (|75p in terms of £(^), we find 



(76) 



(77) 



av 



''coh 



(78) 



The boundary condition is that due to formally infinite 
continuum optical depth if we integrate to v = +oo, we 
have £+ = (where again £ + is the value of £ on the blue 
side of the line). This may fail if the optical depth be- 
tween lines tll is not much larger than 1; the correction 
for this is discussed at the end of Sec. IV CI 

We will concern ourselves with solving Eq. (|78|) in 
Sec. IV Dl using a Monte Carlo method. However before 
we do this we will investigate the implications of the so- 
lution by relating the value of averaged across the 
line, 



f= / t{v)<t>{v)dv, 



(79) 



to the value of A/"^ and to the net decay rate xi\ u , which 
is the quantity we need to know in order to solve helium 
recombination. 



C. Relation of A/^ and £ to decay rate 

The objective of this section is to determine the rate 
±j| u in terms of quantities such as xi, x u , and /; that can 
be computed easily in the level code, and the quantity £ 
that emerges from the solution to Eq. (f78|) . 

As noted at the end of Sec. IV Al the upward and 
downward transition rates in the / <-> u line are Ri u = 



{g u / gi)A u _iJ\f and R u i = A u _>i respectively. Averaging 
Eq. (|77|) over the line profile, we find that 

R lu = —A u ->i [Ac + (A^ 0) - A/c)£]. (80) 
9i 

It follows from Eq. ([65]) that 



+™ Xl A u ^[M c + (M[ a) - A/c)£] }■ (81) 
9i 



The term u x iRiu also appears in Eq. (|7"4")) . so we 

can replace it with M^ . Noting also that A u ^i/T u = 
h = /coh, we may write 

x u = ^ f/i„cA^ 0) + f coh [M c + (M[ 0) -Mc)S]\ ■ 
9i 1 J 

(82) 

The quantity in braces is Al = gix u /g u xi, i.e. the phase 
space density of photons that would be in equilibrium 
with the I and u level populations (again assuming x u <C 
xi): 

Al = /i„cA^ 0) + / coh [Ac + (A^ 0) - Kf c )t] ■ (83) 

If however xi and x u are known (they are directly avail- 
able in the level code), we may rearrange this equation 

to solve for M^ in terms of A/l, Ac, /coh, and £. After 
some algebra, we find 



MP = 



M L - / coh (l - QMc 

1 - /coh(l - 



(84) 



In order to determine the upward transition rate, we 
need to calculate M, which is 



M = Mc + (A/£ 0) - M c )t = M C + 



Then we have 



(A/l - Mc)£ 

l-/coh(l-0 



(85) 



xi\ u = -RiuXi + RuiXu = A u ^i ( -—xiM + x, 

\ 9i 

= ^A u ^ l x l (-M + M h ). 
9i 



(86) 



One can eliminate M using Eq. (|85[) . and algebraic sim- 
plification then yields 

xi\ u = —A u ^ixi(Ml - Mc) (l - - . ^ zv 

91 V 1 - /coh(l - 

(87) 

We now recall that Mc is very nearly equal to the black- 
body function since hydrogen is in Saha equilibrium dur- 
ing He I recombination. 



Xi\ u = A u ^t ( X u - —XlMc ) Pet 

91 



(88) 
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where 



1 



i 



(1 - Oh 



1 - /coh(l - 1 - /coh(l - 



(89) 



Before continuing, wc wish to generalize this equation 
to remove one assumption we have made about contin- 
uum opacity. Recall that we have assumed a finite con- 
tinuum optical depth per unit frequency, which when in- 
tegrated to v — > +oo gives formally infinite optical depth 
and hence drives A/+ to its equilibrium value. In prac- 
tice there are times early in helium recombination where 
the continuum opacity within a line is unimportant, and 
the continuum optical depth between the line of interest 
and the next higher- frequency line is not 3> 1. In these 
cases A/+ may differ from A/pi due to feedback. We may 
trivially correct for this by noting that if feedback is im- 
portant then the line-line optical depth tll is of order 1 
or less. In this case, then the continuum opacity within 
an individual He I line is negligible, so in Eq. 1)75 p wc 
have Tj c f» 0. In this case, we are free to replace A/"c 
in Eq. ([75)) with Af+ without consequence. Making this 
replacement in the subsequent equations, in particular 
Eq. ([75]) . we find that £ on the blue side of the line is 
zero since J\f = A/+ there. Hence the boundary condition 
= that we were using earlier applies, and Eq. (|88|) 
is also valid, except that we need to replace A/c — > A/+ ■ 



%l\u — An— tl 



— XlM+ ) P CS c, 

9i 



(90) 



This equation is very similar to the Sobolev rate, Eq. ([5]). 
The only difference (aside from the factor of 1 + A/+ , 
which is irrelevant since A/+ <C 1 for these lines) is the 
Sobolev escape probability has been replaced with P esc , 
the escape probability including partial redistribution. 



D. Monte Carlo method: theory 

So far in this section, we have constructed an equation 
of radiative transfer (Eq. [78]) for the He I n 1 ?" - 1 ! S 
lines, and related the effective line escape probability to 
its solution (Eq. 159")) . There is only one major step left: 
to solve Eq. (|78p for the full partial redistribution plus 
continuous opacity problem. 

A variety of approaches have been taken in the litera- 
ture for solving line radiative transfer equations including 
coherent scattering terms. One approach is the diffusive, 
Fokkcr-Planck approximation j64[ which replaces the re- 
distribution integral (Eq. 170)) with a second-order differ- 
ential operator. This results in a second-order ODE in- 
stead of intcgro-diffcrential equation, which is a substan- 
tial improvement for most numerical techniques. The 
other possibilities are conversion of the equation of ra- 
diative transfer into a linear al gebra p roblem or solution 
through Monte Carlo methods [66|, M, H, H [zO, O, [zl . 
The latter two have the advantage of being usable in the 
Doppler core of the lines, which we expect to be impor- 
tant since for e.g. He I 2 1 P° — l 1 ^ lines, the width of 



the line that is optically thick to incoherent processes 
Az^jne is only ~ 30A^d during most of He I recombi- 
nation. Therefore we have not used the Fokker-Planck 
approach, which we believe is better suited for studying 
the far damping wings of very optically thick lines such 
as H I Lya (The Fokkcr-Planck operator assumes that 
many scattering events transport a photon over a region 
where the line shape varies slowly; yet in the core of the 
line, single scatterings can transport a photon over the 
width of the core.) We have chosen the Monte Carlo ap- 
proach here, mainly because we had a pre-existing code 
that was capable of handling the problem with minor 
modifications [66| . 

The basic plan of the Monte Carlo simulation is as 
follows: we begin by injecting a photon with frequency 
distribution drawn from the Voigt line profile, (j){v). Wc 
simulate its fate by assuming that it redshifts at the rate 
v = —Hv u i] undergoes coherent scattering with probabil- 
ity per unit time Hv u iTsfi<p(v); undergoes continuum ab- 
sorption with probability per unit time Hv u irj c \ and un- 
dergoes incoherent absorption with probability per unit 
time Hv u its(1 — fi)<p{y). The simulation is terminated 
if the photon is absorbed in the H I continuum, if it 
redshifts out of the line, or if it undergoes incoherent 
absorption in He I. Note that within the idealized con- 
ditions r) c = constant of the simulation, a photon that 
redshifts out of the line will eventually be absorbed by 
the continuum so long as rj c > (though in reality the 
photon would eventually reach other He I lines or redshift 
to below 13.6 eV). Therefore only the total probability of 
these two results is meaningful. We thus denote by Pmc 
the probability that a photon in the Monte Carlo is ter- 
minated by redshifting out of the line or by continuum 
absorption, and let 1 — Pmc denote the probability that 
the photon is terminated by incoherent absorption. The 
implementation of the Monte Carlo is described in Ap- 
pendix (Cj the rest of this section will be devoted to the 
problem of extracting £ from the Monte Carlo. In par- 
ticular, wc will prove that £ = 1 — Pmc- 

Our proof goes as follows. We begin by considering the 
probability distribution n(i/) of the photon frequency in 
the Monte Carlo at any specified time. [Note that this is 
not the same as the histogram of frequencies at which the 
photon scatters, which is cx H(v)(f>(i/) because the scat- 
tering rate is proportional to <p{y) .] Now from elementary 
considerations the probability distribution satisfies 



tl(u) = Hv ul \ — - J7 c n(i/) - T inc( f){v)Tl{v) 



7~coh 



i/)n(i/)- / (t>(v')n(v')p(v\v ') dv 



-T M 4>(v), 



(91) 



where T-^ is the rate of injection of photons in the Monte 
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Carlo. Then the steady-state distribution satisfies 



Now recalling Eqs. (|94|) and ([95]) . we see that 



dll 

— = rjjl{v) + T mc (p{v)IL(u) 



7"coh 



_Lmj_ 

Hv„i 



(92) 



This equation has a striking resemblance to Eq. (|75|) . 
They are both linear inhomogencous integro-differential 
equations, and they obey the same boundary condition, 
namely that £ and II go to zero at high frequency (the 
photons are in equilibrium with H I on the blue side of 
the line). The only difference is in the source terms: the 
source term in Eq. (|78|) is — T mc (j){v), whereas the source 
term in Eq. (|9"2")) is — (T^/ Hv u i)(f)(i/). These source terms 
differ only by a constant scaling factor, and hence the 
solutions are scaled versions of each other with the same 
scaling factor: 



T-, 



HV U IT U 



(93) 



We now note that in the Monte Carlo, photons are in- 
jected (the simulation is restarted with a new photon) 
when the photon is absorbed in cither an incoherent pro- 
cess or by H I. Practically, the frequency span of the 
simulation is finite, yet photons that eventually rcdshift 
through the simulation boundaries will eventually be ab- 
sorbed in a continuum event (assuming the opacity is 
non-zero). Therefore Ti n j = rj nc + r cont , where r; nc 
and r cont are the rates of removal of photons by inco- 
herent scattering and continuum absorption respectively. 
Of these, rj nc is obtained by averaging the incoherent 
scattering rate Ti nc (j)(v) over the probability distribution, 
which yields 



Tine = Hv u iTi nc Il. 



(94) 



[Here ft denotes averaging of Ii{v) over the line profile 
4>{v) , which is equal to the average of the line profile over 
the photon probability distribution II(z/).] The rate of 
removal via continuum opacity is simply 



Tcont = Hv u ir) c . 
Therefore we may write Eq. (|93|) as 

n(^) = (n 



(95) 



'/c 



(96) 



Multiplying both sides by <j){v) and integrating yields 



n = n 



from which we may find 



n 



n + i] c /tu 



(97) 



(98) 



rv 



I-Pmc (99) 



This result connects the Monte Carlo algorithm to the 
parameter £ needed in the level code. Written directly in 
terms of Pmc , the modified escape probability (Eq. [50]) 



is 



P« 



Pmc fit 



1 — /cohP 



MC 



(100) 



Before continuing, we note several properties and lim- 
iting cases of Eq. (|100[) . Since it is obvious that Pmc and 
/coh, being probabilities, are in the range from to 1, 
and we recall that / co h + /mc = 1, it is easily seen that 
P esc is also in the range from to 1 (hence its interpre- 
tation as an "escape probability"). In the limit where 
there is no coherent scattering, we have, unsurprisingly, 
Pesc = Pmc: the effective escape probability P csc appear- 
ing in the rate equations is exactly the probability of line 
escape or continuum absorption in the Monte Carlo. In 
an optically thin line where essentially all photons in the 
Monte Carlo escape, Pmc = li then we find P csc = 1. 
Finally, if Pmc <C 1 so that almost all photons emit- 
ted in the line undergo incoherent re-absorption, we have 
Pcsc = Pv[c/inc : the escape probability in the rate equa- 
tions becomes the escape probability in the Monte Carlo 
times the fraction of the Sobolev optical depth due to 
incoherent processes. 

The Monte Carlo method described here raises two the- 
oretical questions: first, does the escape probability in 
Eq. (|100j) coincide with that calculated in Sec. lIVl as one 
would expect; and second, is the standard Sobolev result 
recovered in the absence of continuum absorption? In 
Appendix [D] we show that the answer to the first ques- 
tion is in the affirmative. As for the second question, we 
find (also in Appendix|D]) that by setting i] c = the stan- 
dard Sobolev escape probability is recovered if the line 
is sufficiently optically thick, i.e. if Pmc "C 1 and if the 
probability P + that a photon entering the line from the 
blue side will escape to the red side without undergoing 
any incoherent scattering satisfies P+ <C 1. 



E. Incorporation in the level code 

In principle, modifications to transport in the entire 
n 1 P° —l 1 S series, the quadrupole series {n}D — 1 1 S'), and 
intercombination series (n 3 P° — l 1 S) of He I could lead to 
acceleration of He I recombination. To make the prob- 
lem computationally tractable, the hierarchy of higher 
order n contributions needs to be cut off where there are 
diminishing returns. The integral solution to the trans- 
port equations with complete redistribution (Sec. IIV[) . 
while not accurate for n 1 P° — l x S rates, gives a quick 
estimate of whether, for example, continuum effects be- 
come more and more important for higher n, or die away. 
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We evaluate the probabilities at z = 1606, at the end 
of He I recombination where the continuum effects are 
expected to be largest, using the integral Eq. (|53[). It 
is found that in the intercombination lines, continuum 
effects become small for n > 4 (the ratio of the proba- 
bilities with and without continuous opacity approaches 
1). For quadrupolc transitions (n 1 !) — l 1 ^), the correc- 
tions are significant out to moderate n, (modifying the 
escape probability by a factor of ~ 2 at n — 8). In the 
— 1 X S series, continuum effects lead to significant 
modifications to the escape probability (by a factor of 
~ 100 at n = 9 toward the end of He I recombination)- 
the entire n 1 P° — l 1 ^ series sees significant corrections. 
(One mitigating factor that makes the calculation prac- 
tical is that only the 2 1 P° — l x 5 line dominates recom- 
bination rates.) Based on these estimates, we calculate 
corrections to transport within the line up to n — 6 in 
the v^P" — \ X S series, up to n = 4 in the intercombi- 
nation lines (n 3 P° — l 1 ^), and n = 6 in the quadrupolc 
lines (n^D — 1 J 5). The n l P° — \ X S series is treated via 
the Monte Carlo technique, while the others are treated 
here using the integral method of Sec. II VI (Fig.llOlshows 
that higher allowed transitions are negligible. In Paper 
III we will show results where the Monte Carlo method 
was used for intercombination and quadrupole lines as 
well.) Fig. ITTI shows the the escape probabilities subject 
to H I opacity (with and without coherent scattering) for 
2 1 P° — l 1 ^ derived from the Monte Carlo, relative to the 
ordinary Sobolcv results. 

Both the Monte Carlo and analytic integral methods 
of finding the escape probability are prohibitively slow to 
run in real-time with the level code. In the Monte Carlo 
method, computing the large number of coherent scatters 
for one photon trajectory in the He I n 1 P° — l 1 ^ series is 
very time-consuming. In the analytic integral method de- 
veloped for complete redistribution, the variety of scales 
in the line profile gives, generally, integrands with large 
dynamic ranges that need to be known to high accuracies. 
For this reason, we generate tables of the escape proba- 
bility over a range of parameters and log-interpolate to 
find the probabilities in the level code. 

The escape probability is a function of the coherent, 
incoherent, and continuum optical depths, and the mat- 
ter temperature, which sets the Gaussian width for the 
line, 

Pesc ( { 7"coh j 7"inc ) ?7c, T m , Tline}) (101) 

The wing width is set by rii nc . We will project to a more 
natural parameter space for the level code, 

Tine? ?7c; ^m? Fling} — > < Z, ,XHeIf, (102) 

I ™HI,Saha J 

since for a given cosmology the latter parameter set de- 
termines the former. We calculate tables for 11 linearly- 
spaced z values between 1400 and 3000, inclusive, and 
21 logarithmically-spaced xjjei values from 2 x 10~ 5 to 
0.08. When we double the fineness of the probability grid 



over the parameters, the change in free electron density 
|Ax e | has a maximum of ~ 5 x 10 -4 at z ~ 1900. The 
neutral hydrogen population is taken to be the evolu- 
tion determined by the reference level code. This is very 
nearly Saha until the end of He I recombination, by which 
time we find that He I has been relaxed to equilibrium. 
(That is, whether you assume the H I population is Saha 
or evolves through in the full recombination treatment 
should matter little for the evolution of the He I ground 
state: even before neutral hydrogen departs significantly 
from equilbrium, He I is already relaxed into equilibrium, 
by that point.) This set of probabilities required ~ 4 days 
to evaluate over 50 x 3.2 GHz computer nodes. The con- 
vergence criterion is described in Appendix[Cl and shown 
to converge to 1.25% fractional error in probability with 
negligible bias. Doubling the number of Monte Carlo 
trials led to a maximum change of |Ax e | < 10~ 4 , and 
resampling with a different random generator (Numer- 
ical Recipes ran2 instead of ran3 |58j |) gives a change 
|Ax e | < 1.5 x 10~ 4 . (More convergence tests are de- 
scribed in Paper III.) 

The result of including these in the level code is shown 
in Fig. [T2J Here, we can test the diminishing-returns 
set of modified lines by running a level code with: 1) 
just 2 1 P° — l 1 ^ modifications, 2) modifications to the 
n 4 P° - 1 X S series up to n = 6, 3) 2 4 P° - 1*5 and the 3 
intercombination lines, and 4) the n 1 P° — l 1 ^ series up 
to n = 6, plus the three intercombination lines, shown in 
Fig. [TUJ It is apparent that most of the effect is due to 
modification to the 2 1 P°-1 1 S" escape probability, leading 
to relaxation of the l n = 2' bottleneck and further that 
coherent scattering is only a small contribution to the 
recombination history. This greatly improves prospects 
for including continuum effects in a practical level code 
through modifications to 2 4 P° — l 1 ^ under the approxi- 
mations developed in Sec. [TV] (without coherent scatter- 
ing). 



VI. CONCLUSION 

The efforts of a new generation of precision small- 
scale CMB temperature anisotropy experiments need to 
be complemented by confidence that the underl yin g the- 
ory is well-understood. Recent results [13, [H, |39T | have 
brought up the possibility that new effects may modify 
the H I and He I recombination histories. Here we have 
examined and extended several of these recently proposed 
effects: the matter temperature evolution and the inter- 
combination and quadrupolc lines of He I. We also intro- 
duce two new effects: feedback of radiation between lines 
(non-local radiative transport), and the effect of contin- 
uum opacity from photoionization of neutral hydrogen on 
transport within and between lines. The effect of feed- 
back between lines, which has previously been neglected 
is also shown to be a significant effect at the ~ 1% level 
by using an iterative method to include the effect. (Pa- 
per II and III will treat two-photon processes, electron 
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FIG. 10: Modification to the He I recombination history due 
to continuum opacity in the n 1 D — l 1 ^, n 3 P° — l 1 ^, and 
n}P° — 1 1 5 lines (and where coherent scattering through 
2 1 P° — l x S is neglected) relative to a model where only con- 
tinuum opacity in the 2 1 P° — PS line is accounted for. The 
total effect on x e is at the level of 1CP 4 . Without feed- 
back, n P° — 1 1 S slightly speeds recombination relative to just 
2 1 P° — 15— the allowed rates from n > 2 to the ground state 
are accelerated. Coherent scattering through 2 1 P° — PS is 
also seen to be negligible for the recombination history. 
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FIG. 11: The modified escape probability from 4 He 2 1 P° - 
PS during He I recombination, comparing the results of the 
standard Sobolev approximation and the modification due 
to continuous opacity with and without coherent scattering. 
Once coherent scattering is introduced, photon diffusion ef- 
fects increase the escape probability by increasing the span 
of frequencies traversed by the scattering photon before it es- 
capes. These probabilities are log-interpolated over the grid 
of XKei and z used in the level code. We find that the effect 
of doubling the grid resolution and with it the smoothness of 
the interpolated probability is negligible. 



and 3 He scattering, and rare decays.) 

The striking new effect has been the introduction of 



FIG. 12: Modifications to the He I recombination history from 
the inclusion of feedback and continuum opacity in He I lines, 
compared to a Saha He I recombination history and the stan- 
dard He I recombination history. Continuous opacity starts to 
become important at z ~ 2100 (as suggested by Fig. [7J), and 
pushes the He I evolution to a Saha by z ~ 1800. The begin- 
ning of H I recombination is visible here starting at z ~ 1700, 
and for later times x e drops precipitously. The modifications 
to He I recombination suggested here are twofold: at early 
times during He I recombination, feedback slows recombina- 
tion, and at later times, continuous opacity accelerates recom- 
bination relative to standard models. 



continuum opacity due to H I photoionization in the He I 
2 1 P° — l 1 S transition. We presented two methods of un- 
derstanding this physics, an integral solution to the ra- 
diative transport equations for complete redistribution, 
and a Monte Carlo method for partial redistribution (to 
include strong coherent scattering effects) . This is found 
to dramatically increase the escape probability and speed 
up He I recombination once continuum effects in the line 
become important. These effects, in combination with 
the inclusion of the intercombination line as suggested 
by paint a very different picture of He I recombi- 
nation - one where the recombination accelerates after 
z ~ 2200 leading to very little overlap with the H I re- 
combination. The speed-up is a generic feature of any 
effects that facilitate the formation of ground state He I, 
though of course the details depend on the kinetics of the 
specific mechanism. 

Now that the interpretation of many important cosmo- 
logical parameters depend on slight shifts (such as evi- 
dence for inflation from the primordial slope n s ), it will 
be crucial to understand even sub-percent corrections to 
the recombination history from various slow processes. 
We have examined several corrections here and more are 
considered in Paper II and Paper III. Once a set of cor- 
rections to hydrogen and helium recombination are so- 
lidified, the next step will be to parameterize a set of 
corrections to a recombination model in a fast, easy-to- 
use package that can be plugged in to CMB calculations 
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and parameter estimations. Distilling the modified es- 
cape probabilities for a set of cosmological parameters 
into a new recombination code will be the subject of later 
work. 
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APPENDIX A: ATOMIC LEVELS AND 
RADIATIVE RATES 

1. Hi model atom 

The full model H I atom has 445 levels up to a maxi- 
mum principal quantum number n max = 400. The indi- 
vidual Z-sublevels are resolved for n < 10. The energies 
of each level are determined from the exact nonrelativis- 
tic formula E n i = —Rn/n 2 , where the hydrogen Rydberg 
constant is Rn/h = 3.2881 PHz. The degeneracy factors 
are g n i = 2(21 + 1) for the /-resolved levels and g n = 2n 2 
for the unresolved levels. Throughout, we have used a 
smaller set of 245 levels, identical to the full model, ex- 
cept truncated at n max = 200. Inclusion of the remaining 
levels up to n max = 400 will be important for late-stage 
H I recombination, but negligible for He I recombination. 

The bound-bound electric dipole Einstein coefficients 
between /-resolved levels are determined by integration 
of the exact wave functions. For cases where the upper 
sublevel is unresolved (n M > 10), we compute a weighted 
average of the Einstein coefficient, with the weights de- 
termined by the degeneracy factors. This gives the ap- 
propriate decay rate in the limit where the sublevels are 
populated according to their statistical ratios. For tran- 
sitions with ni > 10 we used the Gaunt factor approx- 
imation of Ref. [73| |. Eq. (1.38), with the correction de- 
scribed in the Appendix of Ref. [74| . Direct summation 
of the Einstein coefficients shows that this yields an error 
of 1.2% for the worst case (n; = 11, n u — ni = 1). We 
have not included the H 1 electric quadrupole lines (Is— 
nd, n > 3) because they are at the same frequency as the 
Lyman lines ls-np and hence do not yield acceleration 
of the H 1 recombination process. (This is not the case 
for He 1.) 

The bound-free cross sections for n < 10 are obtained 
from TOPbasc (75|; for n > 10, the /-sublevels are unre- 
solved and we have used the Gaunt factor approximation. 

The H 1 atom posesses a metastable 2s level, which de- 
cays by emission of two electric dipole photons. We have 



used the two-photon frequency distribution and the spon- 
taneous lifetime estimate Ahi = 8.2249 s" 1 from Ref. [7a ]. 



2. He 1 model atom 

The He 1 model atom has 289 levels up to a maximum 
principal quantum number n max = 100. For all values of 
n we resolve singlet ( "para- He 1" ) versus triplet ( "ortho- 
He 1") levels. The /-sublevels are resolved for n < 10. We 
include only configurations of the form Is nl as the dou- 
bly excited configurations are all unbound. The energy 
levels for the /-resolved levels in He 1 are taken from the 
NIST Atomic Spectral Database, based on Ref. [77]|, ex- 
cept for n = 9,10,Z>6, for which hydrogenic values were 
used. For the n > 10 levels in which the value of / is un- 
resolved, we have used the Rydberg formula (which pro- 
duces as accurate a value as can be considered meaningful 
given that the various / levels are not strictly degenerate). 
The degeneracy factors are g n LS = (2£ + 1)(25 + 1) for 
L-resolved levels and g n s = n 2 (2S+l) for unresolved lev- 
els. Note that in ortho-He I, there are multiple allowed 
values of J for given L and S, which correspond to fine 
structure levels; we do not resolve these. 

For the allowed bound-bound transition rates, we 
combined data for several sources. For transitions be- 
tween two /-resolved levels, the electric dipole oscillator 
strengths for allowed S — P° and P° — D transitions were 
obtained from Ref. [78[ for upper levels with n < 9. The 
l 1 ^ — lO 1 /-" 3 transition rate is from Ref. [79| and other 
S — P and P — D rates with upper level n = 10 are 
from TOPbase [z|. For D - F°, F° - G, etc. transi- 
tions, we used the hydrogenic rates since the I > 2 states 
are well approximated as an electron orbiting a pointlike 
He + ion. For transitions between an Z-resolved lower level 
and an unresolved upper level (n > 10), we first compute 
the fully /-resolved Einstein coefficients and then average. 
These /-resolved rates are taken as hydrogenic for P — D, 
D — F, etc. transitions. For l 1 ^ — n l P° transitions 
(n < 10), we use the asymptotic formula of Ref. [80| . 
(This formula is only 1.4% different from TOPbase at 
n = 10.) The Coulomb approximation [8l| is used for 
all other S — P° transitions. For transitions between two 
unresolved levels (n; > 10), we have used the hydrogenic 
Gaunt factor approximation. 

In He 1 recombination, the intercombination and for- 
bidden lines can be competitive with allowed lines as 
a mechanism of decay to the ground state because the 
allowed lines become very optically thick. In particu- 
lar we consider the series of intercombination lines He 1] 
l^-n 3 ^ (n > 2) and forbidden lines [He 1] ^S-r^D 
(n > 2). For the He 1] l 1 ^ - n 3 P? lines with n = 2, 3 we 
have used the Einstein coefficients estimated by Ref. [82| ; 
these are extrapolated as A cx n~ 3 for larger n. (Note 
that the Einstein coefficients are usually quoted for the 
fine structure-resolved n 3 Pf level. Since our code does 
not resolve the fine structure these coefficients must be 
divided by 3 since only 3 of the 9 n 3 P° states have J = 1.) 
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For the [He i] l 1 ^ — r^D lines with 3 < n < 6, we have 
used the oscillator strengths calculated by Ref. [&|- For 
n > 6 the oscillator strengths have been extrapolated 
with the asymptotic expansion / oc n~ 3 . 

Like H I, the He i atom posesses a metastable sin- 
glet state 2 1 £o that decays by emission of two elec- 
tric dipole photons. The two-photon rate used here is 
A Hc i = 50.94s -1 , as determined by Ref. Q. The two- 
photon spectrum of He I was determined by the following 
formula, which was fitted to the results of Ref. 1841 : 



19.602 C 3 (l-742 - 7.2C+12.8C 2 ) 



^l 1 s -2i s 



(C + 0.03) 5 



(Al) 



where ^iis _2 1 So = 4.9849 PHz is the frequency differ- 
ence between the ground and metastable states, and 
C = K^iiS -2iS - v )/ v vs -2is is between and 1/4. 



Monte Carlo methods developed in the next section for 
the n 1 P° — l 1 S series, before the effect of coherent scat- 
tering is "turned on" . This requires solving Eq. (|53p with 
the Voigt profile. (It also shows directly why complete re- 
distribution is a dangerous assumption in the n l P° — 1 J 5 
series.) Third, He I recombination rates are sensitive to 
the probabilities set by Eq. (fS"3"|) . so a high-precision nu- 
merical result is necessary. These constraints are met by 
directly integrating Eq. (|53|) with the Voigt profile. 

The Voigt profile presents two characteristic frequency 
scales: slowly varying functions in wings, and rapidly- 
varying functions in the core. In the Doppler core, we 
can follow an analysis similar to [63| and note that in the 
inner integrand, 



exp 



4>(y)dy 



(B2) 



3. He II model ion 

The He n ion has one electron and hence in nonrel- 
ativistic theory all rates can be obtained by scaling of 
H i. The energy levels are re-scaled in proportion to 
Z 2 /i, where Z is the atomic number and [i is the reduced 
mass. The ratio of these factors for He n versus H I 
is 4.001787. The bound-bound rate coefficients scale as 
(Z 2 fi) 2 and the photoionization cross sections scale as 
(Z 2 /j,) -1 . The two-photon decay rate A scales as (Z 2 fi) 3 . 
The He n model atom is Z-rcsolvcd up to n = 10 and has 
a maximum principal quantum number n max = 100, for 
a total of 145 levels. 



APPENDIX B: INTEGRATION OF THE 
TRANSPORT EQUATION FOR COMPLETE 
REDISTRIBUTION AND CONTINUOUS 
OPACITY 

In this Appendix we discuss how the transport equa- 
tion with continuous opacity is integrated to find the con- 
tinuum absorption probability, Eq. (|53[) . Rybicki and 
Hummer [63| develop a convenient analytic approxima- 
tion to the transport solution in the case of a Gaussian 
profile with high optical depth. This gives good intu- 
ition for the combinations of factors that contribute to 
the escape probability. For example, it suggests that if 
the escape probability for a line i is known, then an ap- 
proximate scaling (in their limits) can be used to find the 
escape probability from j , 



-Pesc,i7inc,M 0"D,j?7C,j 



OD,i?7C,i 



(Bl) 



'mcj 



There are three reasons why we would like a more general 
solution. First, the high optical depth limit is not ap- 
plicable to intercombination or quadrupole lines during 
recombination, though they have nearly Gaussian pro- 
files. Second, it is desirable to be able to use the macro- 
scopic transport equations as a cross-check for the photon 



the contribution is small unless the integral in the ex- 
ponent is less than t^. Yet, for the small differences 
between v and v, the integrand is nearly linear, 



(j>{v')dv' w (j){v)(i> - v). 



(B3) 



Thus, in the Doppler core, the integral depends on 
evaluations of the integral over the Voigt profile on very 
small scales. The Voigt function and its integral are time- 
consuming to evaluate. Further, a linear interpolation 
lookup table must have sufficient resolution in the core, 
and breadth in the wings for Eq. (|53p to be evaluated 
accurately. We use the Gubner's series [85| in the core 
region and switch to a fourth order asymptotic expan- 
sion of the Voigt profile in the wings. (In the bound- 
ary between the two regimes, we estimate the error in 
the asymptotic expansion by using the next higher order 
and ensure that differences between the two approaches 
are negligible.) The Voigt function and its integral are 
tabulated out to 10000 Doppler widths with 50 values 
per doppler width, exploiting the symmetry of the pro- 
file and its integral. These are interpolated using a cubic 
spline and Eq. (f5"3")) is integrated using a 61-point Gauss- 
Kronrod adaptive integration scheme. (This is necessary 
because the integrand of Eq. (f53|) carries its support over 
a widely varying range of v.) To apply the numerical 
results of this integral to the recombination setting, Pq 
is pre-calculatcd for a range of r; nc and r] c Ai/D^ and then 
log- interpolated. 



APPENDIX C: MONTE CARLO METHOD: 
IMPLEMENTATION 

Monte Carlo radiative transfer methods have been used 
extensively in the literature [H, 13, H, [110, [lH, [z3| - 
We follow an approach similar to Hirata 66l|. and use 
an algorithm to draw atomic velocity distributions that 
is discussed clearly in Lee J8f| H?}- A general atomic 
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scattering process can be represented by the joint prob- 
ability function of a photon with v m scattering against 
the atom coherently to produce a photon of u ont through 
the scattering angle x, P(v<mt,xWin)- 

In coherent scattering, the photon is emitted with the 
same frequency it was absorbed in the atom's rest frame, 
and its emission direction relative to the incoming photon 
is drawn from the angular distribution associated with 
the transition. The Doppler shift formula gives an ex- 
pression for the frequency shift between incoming and 
outgoing photons in the lab frame, 



= (/|| + a)(l - cos x) ~ f-L sin x, 



(CI) 



where f = z/ v/c, for the transition frequency is , and 
a = hv" 2 / {m p c 2 ) . Thoughout, || labels the component of 
a vector in the direction from which the photon came, 
and _L labels the component perpendicular to this direc- 
tion and in the plane of scattering. Thus, the change in 
frequency between the input and output states in the lab 
frame is uniquely specified through Xi f\\ an d /_!• These 
three quantities are stochastic, where f\\ and f± depend 
on the thermodynamics of the gas, and x is determined 
by the quantum mechanical scattering distribution. The 
distribution of the angle x between outgoing states and 
incoming states for dipole scattering is given by (over the 
range < x < [66j . 



p (x)dx = o ( X + cos2 X) sin X dx- 



(C2) 



We can scale the atomic velocity relative to the other 
velocity scale in the problem, the characteristic thermal 
velocity, u = v y/m/(2k B T). This gives the convenient 
expression = (Aisd)u\\ = ^/2<tdU\\, where A^d is the 
Doppler width in standard notation, and a 2 D is the vari- 
ance of the Doppler Gaussian. In the perpendicular di- 
rection, the distribution of atomic velocities is thermal, 

2 

oc e~ u± . Bayes' rule gives the distribution of ttii based 
on known distributions as 



P(u\\\x in ) 



P(x in \u\\)P(u\\) 
a 



(C3) 



nH(a,x in ) a 2 + (x m - uii) 5 



Here H(a, x) is the Voigt profile, 



H(a, x) 



dy 



(C4) 



x = (y — vq)/Ai/£>, and a is the Voigt width parameter 

Inline 



47rA^ D ' 



(C5) 



Typically n«l, indicating that the Lorentzian width of 
the line is very small compared to the Doppler width. Lee 
[HI, [8^ | presents a rejection method to randomly draw 
from this distribution. 



When a photon is emitted by an incoherent process - 
i.e. it reaches an excited level u either by a radiative tran- 
sition from another excited level, or by recombination - 
the outgoing photon has no knowledge of the phase space 
distribution of pre-existing line photons, so the probabil- 
ity distribution function for an outgoing photon is just 
the Voigt profile. Because the Voigt profile is the con- 
volution of the Doppler Gaussian distribution and the 
Cauchy distribution, a Voigt-distributed random num- 
ber is easily implemented as the sum of random numbers 
drawn from each distribution. 

The distance between scatters for small frequency 
shifts is given by I f=a ci/ _1 (z)^ 1 Av (for some central 
frequency v$ and shift A^). Over the travel times asso- 
ciated with the escape of one photon, the universe has 
expanded very little, and the line depth parameters are 
effectively constant. 

In standard Sobolev theory, where the fate of the pho- 
ton is either that it escapes or is absorbed in an inco- 
herent process, it is clear what is meant by an escape 
probability. The transport problem for He I with con- 
tinuum opacity is not as clear-cut. Coherent, incoherent, 
continuum and redshift processes are all active, so there 
are several choices about what "escape" and "scatter" 
mean. As shown in the main text however, there is one 
specific number we need to know: given a photon was just 
emitted in the line through an incoherent process, what 
is the probability Pmc that it will escape through red- 
shifting or continuum absorption before it is absorbed in 
an incoherent process? The easiest way to measure Pmc 
with the Monte Carlo is to inject many photons whose 
initial frequency distribution is the Voigt profile, and fol- 
low them until they escape (i.e. redshift out of the line 
or get absorbed by H i) or are re-absorbed in an incoher- 
ent process. Here the range of frequencies simulated is 
±360 THz with a bin size of 0.36 GHz. (Using a smaller 
span or less resolution leads to escape probabilities that 
are biased high at early times during He I recombina- 
tion, because absorption can occur far in the wings at 
early times when the incoherent depth is high and the 
continuum depth is small. Doubling the boundary and 
halving the frequency step size does not improve the re- 
sults, within these tolerances.) The basic steps in the 
Monte Carlo arc: 

1. Draw a photon from the Voigt distribution, repre- 
senting emission from an incoherent process. 

2. Draw the optical depth that the photon traverses 
before scatting from the exponential distribution, 
P{5t)(1{5t) = e- Sr d(5T). The next frequency 



where the photon scatters is given implicitly by 

(7"inc 



dv 



dvr)(v), 



(C6) 



where we have identified the integrand dr tot / dv as 
I(v). This is implemented numerically by choosing 



25 



frequency bins (Au = 0.36 GHz) and determining 
whether a photon crosses that bin, or is absorbed. 
If the bins are chosen to be small enough, then the 
integrand is nearly linear over the bin, 



Vi+i - Vi 
Ay 



Vi- 



(C7) 



Integrating this gives a quadratic expression for the 
fraction x of the bin that the photon traversed: 



Vi+i 



2rn 



Vi 2 
X 



St 



{Av)rj t 



(C8) 



This can either be solved quadratically, or by work- 
ing to first order in small bins. We will use the 
latter approach for small bins, 



x 



St 



{Av)r}i 



and x = xq ( 1 — xq 



Vi +i ~ Vi 

2m 



(C9) 



The fraction of the bin traversed indicates whether 
the photon scattered in the bin, or leaves. If it es- 
caped the bin, then we move the photon to the start 
of the next bin and goes to step 2. If it traversed 
less than the whole bin, we find the frequency in 
the bin where it scatters v S catter, and g° to step 3. 

3. Draw a uniform random number between zero and 

(^scatter) and determine the type of event (inco- 
herent scatter/absorption, coherent scatter, or H I 
continuum absorption). 

4. If the photon coherently scatters, draw the scat- 
tered atom's velocity and angle between incoming 
and out-going states and use Eq. (|C1[) to find the 
photon energy after scattering, and go to step 2. 
If the photon is incoherently scattered by He I, un- 
dergoes photoelectric absorption by H I, or redshifs 
through a pre-defined simulation boundary, start a 
new photon in step 1. 

The Monte Carlo procedure leaves one issue open, 
namely the convergence criterion. We repeat the Monte 
Carlo until 6400 photons escape from the line, which 
should enable Pmc to be determined to a fractional error 
of6400~ 1 / 2 = 0.0125. A possible concern with this proce- 
dure (or any other in which the number of photons simu- 
lated is not fixed before running the Monte Carlo) is that 
the result could be biased if the convergence criterion de- 
pends on the results of the simulation. We addressed this 
question by replacing the function that decides whether 
the photon escapes with a random number generator that 
returns "escape" a known fraction Pmc of the time and 
"no escape" (incoherent scatter) a fraction 1 — Pmc of the 
time. This produces a very fast code that allows us to 
map the distribution of the estimated escape probability 
Pmc as a function of the true Pmc- Across the relevant 
range of Pmc (down to 10~ 6 ) we find that Pmc has a bias 
(Pmc)/Pvic — 1 whose absolute value is < 0.1%, and a 
standard deviation ct(Pmc)/Pmc ~ 0.0125. 



APPENDIX D: RELATION OF ESCAPE 
PROBABILITY METHODS 

In this Appendix, we consider the relation of the Monte 
Carlo method to the analytic solution of Sec. IIVI and to 
the usual Sobolcv escape method. We show explicitly 
that the solution to the transport equations considered 
here reduce to those of Sec. IIVI if / co h — ► 0. We then 
proceed to investigate the conditions under which the 
Monte Carlo solution is equivalent to the Sobolev escape 
method. 



1. Relation to case of no coherent scattering 



Here we connect the formalism developed in Sees. IV Bl 
and IV CI to the analysis of Sec. IIVI which neglected co- 
herent scattering. 

The transport equation for line processes and the con- 
tinuum, ignoring the coherent scattering diffusion term 
gives 

^ = Vc [M{v) - Afc] + n nc j>{v)[N{v) - Af[ 0) }. (Dl) 

This can be solved for tf analogously to the derivation 
in Sec. HVl 



(D2) 



where Pesc = Ps — AIl is the escape probability de- 
rived in Sec. IIVI (c.f. Eq. l55|) . Comparing this to the 
analysis of Sec. IV C| in particular Eq. (|76|) we find the 
correspondancc 



1=1- P( IV ) 

S x A esc ' 



(D3) 



Plugging in to the transport result, Eq. one finds 

that P esc (which modulates the rate) is now 



p — 



p(IV) f 
-'esc J i 



1 - f ,P (IV) 
1 Jcoh^csc 



(D4) 



This reduces to P csc = Pelc in the limit that scattering 
in the line is purely incoherent, as expected. 



2. Relation of Monte Carlo to Sobolev method 

Next we consider the relation between this escape 
probability and the traditional Sobolev probability Ps 
used in recombination codes in the absence of coherent 
scattering. The Sobolcv escape calculation in its usual 
form assumes complete redistribution over the steady- 
state line profile and neglects continuum opacity [64], l8a | . 
Complete redistribution is an accurate assumption in 
pure-resonance rate problems where the radiation field is 
forced into equilibrium with the line over most of its cen- 
tral frequency extent by the high incoherent scattering 
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rate (typically when the wing is optically thick to inco- 
herent scattering). Also it is trivially valid if the scatter- 
ing optical depth is dominated by incoherent scattering. 
Diffusion from Doppler shifts in repeated coherent scat- 
tering events has been studied extensively through its re- 
distribution function, the Fokker-Planck approximation 
(accurate away from the Doppler core), and Monte Carlo 
methods [6J, [66j, [86j, |88| . This diffusion tends to broaden 
the jump in the radiation phase space density distribu- 
tion on the blue side of the line, but does not significantly 
suppress (or enhance) Af near line center or modify the 
atomic transitions rates. Thus, for practical purposes in 
a recombination level code with pure resonant line pro- 
cesses and no continuous opacity, there is very little loss 
in assuming complete redistribution and steady-state ra- 
diation fields. 

During recombination, there are two general categories 
of lines connecting to the ground state to consider: (1) 
fine ~ 1 so that Ti, lc Ri rg (intercombination and for- 
bidden lines) or (2) /i nc is small but r; nc is large, i.e. 
Tcoh ^ Tine ^ 1 (allowed lines such as n}P° — IS in 
He I or the Lyman series in H i). We consider what 
happens in each type of line with continuum opacity off 
(Vc = 0). 

For case (f), / co h <C 1 but /i nc « f; then it is imme- 
diate that P esc ~ Ps- The treatment of case (2) is more 
complicated. We begin by integrating Eq. (|75|) over fre- 
quency with rj c = to get 



£+ — £- = T coh 



4>{v'W)p{v\v')dv' dv 
4>{y)\i(y)-\\dv. 



(D5) 



Here £± is the value of £ on the blue (+) or red (— ) 
side of the line. The term multiplying T coa vanishes since 
the redistribution probability p integrates to unity. Also 
£ + = because of boundary conditions. The last term 
simplifies to Ti nc (l — or Ti nc py[c- If we turn off contin- 
uum opacity (rj c = 0), this then simplifies to 



Pmc = — ■, (D6) 

7"iiic 



which is much less than 1. Then from Eq. (|100|) 
p 

1 CS' 



(£— /uncO/inc ^ /inc (D7) 

1 — /coh-RvTC Tine TS 



In order to proceed we must understand the behavior 
of in the case with no continuum opacity. This can 
be done by considering a thought experiment in which 
we inject photons into the Monte Carlo on the far-blue 
side of the line instead of using the Voigt profile. In this 
modified Monte Carlo, the photon probability distribu- 



tion n mod satisfies 

3n mod 



di 



7c oh 



v)u mod {v) 
<^)n m °» 



v')u mod {v') P (v\v')dv' 



(D8) 



with the boundary condition n mod (;/) = Ti n j/Hv u i in the 
blue wing since in the absence of the line photons rcdshift 
at a rate v = —Hv u i. Defining the quantity 



X{v) = I - 



-n D 



we see that X(-\-oo) = and that X(v) satisfies 



^"coh 



(D9) 



{v)X{y) - J <t>(v')X{v') V (v\v')dv' 
{v)\X{v) - I], (DIO) 



where we have used the symmetry relation 
/ <j)(v')p(v\v') dv' = 4>{v) to eliminate the terms 
that come from the 1 term in Eq. (|D9[) . Now X(v) 
satisfies the same integro-differential equation as £,{v) 
in Eq. ([TS)) . and has the same boundary condition 
X(+oo) = 0. Therefore we have X{v) = £{v) and 
£_ = X_. This means that the photon frequency 
distribution on the red side of the line is 



n 



mod 



■(i-e-) 



(Dll) 



Therefore on the red side of the line, the frequency distri- 
bution of the photons is suppressed by a factor of 1 — £_ . 
Therefore, in the absence of continuum opacity, the prob- 
ability Ptrans that a photon injected on the blue side of 
the line manages to redshift through the line without any 
incoherent scattering/ absorption is 1 — Conversely, 
= 1 — -Ptrans is the probability that such a photon is in- 
coherently scattered/absorbed in the line. This provides 
a simple physical interpretation of . 

In the case where the line is very optically thick into 
the damping wings, we expect the line transmission prob- 
ability Ptrans <C 1 and hence w 1. It follows that 
-Pcsc ~ T s 1 - This is in agreement with the Sobolev value 
for rs ^> 1. Thus, in both of the limiting cases (1) and 
(2) where the continuum opacity is small, P csc is well- 
approximated by Pg. The approximation P csc ps Pg can 
fail if the coherent depth is very large, while the incoher- 
ent depth is optically thin. This could occur at very early 
times during recombination in the n 1 P° — l 1 S series, for 
example. Yet, for 2 1 P° — l 1 ^, the difference between 
/incPnc and P s is < 1CT 4 (likewise for the n x P° - \ X S 
series considered here). The distinction is thus negligible 
in the overall recombination history. 
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APPENDIX E: PHOTONS IN THE HE I 
CONTINUUM 

Recombination to the ground state is usually ignored 
because the photon that is generated immediately rc- 
ionizes another bound atom. This is a single-species 
picture. When bound-free rates in helium and hydro- 
gen interact, the recombination photons that would have 
ionized He I can instead be absorbed by H I. This ac- 
celerates He I recombination by allowing recombinations 
directly to l 1 ^. (It also forces hydrogen slightly out of 
equilibrium.) The direct l 1 ^ recombination rate in He I 
should roughly scale as the recombination rate to l 1 ^ 
times the probability that the ionizing photon is absorbed 
by hydrogen instead of by He I. In this appendix, we de- 
velop this effect in more detail and show that it can be 
neglected. 

It is tedious to solve the full transport problem for 
Af(v) subject to bound-free processes in He I and H I, 
and then integrate to find the photoionization and recom- 
bination rates. We will instead track the total number 
of photons with energies above 24.6 cV (from recombi- 
nation to He i). and develop a method to calculate an 
effective cross section and photoionization rate for both 
He I and H I by photons from region. Recombination 
rates can be calculated through detailed balancing of the 
photoionzation rate. 

A similar mechanism is active in He II, and in princi- 
ple one should also consider simulating the possibility of 
direct recombination to H I Is. In H I, there is no cou- 
pling between species and the answer is straightforward. 
In He II, recombination proceeds too close to equilibrium 
for this effect to matter. 

We lump the radiation field into frequency regions for 
photons above threshold for each species and consider the 
loss of photons in a given lumped region due to bound- 
free processes, and redshift. Note that the integral of 
the radiation phase space density over a region is not 
enough to entirely specify the behavior. The number of 
photons that redshift out of the region depends on the 
phase space density at its red boundary. Because the 
spectrum is black-body dominated, we will approximate 
the radiation phase space density at the boundary by its 
black-body value there. 

If the perturbations to the radiation field caused by 
bound-free processes are over frequency scales which are 
large compared to the exponential fall-off of the black- 
body radiation field and the radiation field is close to 
equilibrium, then an effective constant cross-section (near 
threshold) in each lumped frequency region provides a 
good first approximation. 

Define the number of photons per hydrogen nucleon 
above the ionization threshold of species s as 



dE 



(El) 



Af(i 



dv. 



(Throughout, we will use the subscript s to denote the 
species. Unless stated otherwise in this appendix, all 
occupation variables, photoionization cross sections, and 
photoionization or recombination rates are taken from 
the ground state of that species.) To a good approxima- 
tion, at these energies stimulated recombinations can be 
neglected because the photons are from well into the tail 
of the blackbody distribution. One can then write the 
photoionization rate from the ground state due to these 
photons as 



0s 



Ccff.i 



,cn H 



h 3 



(E2) 



where er ff. s is an effective cross section. By the principle 
of detailed balance the recombination rate directly to the 
ground state is 



Coff,, 



n e n c 



LTE 



X 



>, bb 
th. 



h 3 



(E3) 



The effective cross sections are determined by equating 
Eq. (|E2[) to Eq. (|4|) for a thermal distribution of photons 
Af = e -> w / k * T , i.e. 



Ocff.s 



! /°° v 2 G{v)e~ hv l k * T dv 



a(u n 



X th, s 
i) H 7— CT {Vn 



0: 



(E4) 



where a 1 denotes the derivative of the cross section with 
respect to frequency. Note that this assumes the dis- 
tribution of radiation above the threshold v m in is dis- 
tributed as Af oc e ~ hlJ / k n T . This is obviously true in 
equilibrium if hv 3> k^T but may be violated in the 
real universe. Physically we do not expect large devi- 
ations from this proportionality if the main source of 
opacity is photoionization, because the photons above 
ionization threshold should develop a phase space distri- 
bution Af = e -i hv -^)l k ^ T t where the photon chemical 
potential /i 7 is determined by the non-Saha behavior of 
the ionization stages, e.g. 



M7 = Mhii + Me - A f m = fc B Tln 



n e n m i/n m 
'n e nHii/n H i)saha 



(E5) 



if the source of opacity is H I. This may not happen 
if the photon sees opacity from two competing species 
(e.g. H I and He i) whose relative deviations from Saha 
equilibrium are different. The most severe case would 
occur if the cross sections vary rapidly with frequency 
such that e.g. Hi opacity dominates at some frequencies 
and He I at others. However the H I and He I cross sec- 
tions are smooth functions of frequency in the regime of 
interest (i.e. far below the He I double-excitation reso- 
nance region) , so we expect that the current "one group" 
approach is adequate for assessing whether continuum 
photons are important. 
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FIG. 13: Recombination directly to He I l 1 /? can occur be- 
cause He I recombination photons can ionize neutral hydro- 
gen, before they ionize He I. This would tend to accelerate 
recombination. Here we show that this effect is leads to a 
maximum fractional deficit of electrons of around 10 -4 . 



Then, the total bound-free rate to the ground state is 

LTE 



dx, 
~~dt 



cn H 



h 3 



Tl e Xc 



n e n c 



X>> bb - T X> 



(E6) 

In this picture, in addition to the state occupations evolv- 
ing, the lumped radiation fields X{t evolve, and can be 



included as additional "states" in the level code. The 
total number of photons in a frequency region is not con- 
served. Thus the evolution equations must account for 
photons that redshift out of the region, as well as those 
that are injected or removed by bound-free processes. 
Taking the derivative gives a loss rate from the region, 



dt 



1 



71h dt 
8tt h 3 H 



y th- 



(E7) 



Here vth is the ionization threshold frequency of the 
species. Assuming an Af oc e - hu / k sT S p ec trum, the rate 
of redshifting of photons is (under the considerations de- 
veloped for He I recombination), 



1 



dt 



bb 



hv 
k B T r 

hv 
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H. 



With this choice of variables and tracked states, we 
can accommodate the interaction between He I and H I 
bound-free rates-the overall effect on the free electron 
occupation fraction is shown in Fig. 1131 
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